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PREFACE 


This  report  is  the  Phase  n  Small  Business  Innovative  Research  (SBIR)  Final  Report 
under  U.S.  Air  Force  Contract  No.  F33615-92-C-2288.  The  goals  of  the  project  were  to 
improve  ramburner  fuel-air  mixing  and  combustion  efficiency  through  the 
innovative  design  of  the  fuel  injection /flameholding  system,  and  to  improve  upon 
current  CFD  techniques  for  designing  combustion  systems  by  providing  a  detailed, 
high  quality  database  for  code  validation.  The  project  consisted  of  five  major  parts: 

1)  a  numerical  investigation  of  the  baseline  Integral  Fuel 
injection /Flameholder  (IFF)  design,  and  advanced  designs; 

2)  isothermal  experimental  tests  of  the  baseline  IFF  design,  and  selected 
advanced  designs; 

3)  numerical  combustion  tests  of  the  baseline  IFF  design,  and  selected 
advanced  designs; 

4)  design  and  fabrication  of  experimental  hardware  for  combustion  tests; 
and 

5)  numerical  reacting  analysis  of  an  existing  benchmark  experiment. 

Mr.  Clifford  E.  Smith  was  the  program  manager  and  Mr.  S.  Alan  Spring  was  the 
principal  investigator  and  project  engineer  for  CFD  Research  Corporation  (CFDRC). 
The  U.S.  Air  Force  technical  contract  monitor  was  initially  Ms.  Lourdes  Maurice, 
and  later  Dr.  Abdi  Nejad.  For  the  experimental  tests,  the  U.S.  Air  Force  project 
engineer  was  Mr.  Charbel  Raffoul  and  Mr.  Kevin  Kirkendall  of  CFDRC  helped  in 
taking  the  measurements.  Valuable  assistance  was  attained  from  Mr.  Milind 
Talpallikar,  Dr.  Andy  Leonard,  and  Dr  Yong  Lai  of  CFDRC.  The  final  manuscript 
was  skillfully  prepared  by  Ms.  Marni  Kent. 
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EXECUTIVE  SUMMARY 


The  combustion  efficiency  of  current  ramburner  technology  is  mixmg-hmited. 
Therefore,  significant  gains  are  possible  if  the  fuel-air  mixing  can  be  increased 
without  a  significant  increase  in  pressure  loss  or  a  decrease  in  flame  stability. 
Several  design  concepts  were  explored  in  this  SBIR  Phase  11  project  that  increase  the 
fuel-air  mixing  by  generating  large  axial  vortices  downstream  of  the  flameholder 
The  basXe  flLeholder  for  this  project  was  the  Integral  Fuel  injector/ 
Flameholder  (IFF)  used  in  previous  studies  by  United  Technologies  Research  Center 
(UTRC).  CFD  solutions  were  completed  on  the  baseline  design  and  28  advanced 
configurations.  Most  designs  did  improve  fuel-air  mixing,  but  also  increased  the 
cold  pressure  loss  by  several  fold.  The  most  promising  designs  used  the 
momentum  of  the  fuel  to  contribute  to  the  axial  thrust,  thereby  reducing  the 

pressure  loss. 

Isothermal  experimental  tests  were  conducted  at  Wright  Laboratory  in  parallel  with 
CFD  analysis.  A  complete  dataset  was  taken  for  one  baseline  model,  while  several 
advanced  models  were  screened  for  pressure  loss.  The  baseline  tests  showed  that  a 
regular  pattern  of  alternating  vortices  were  shed  from  the  base  of  the  flameholder. 
Simultaneous  three-component  laser  Doppler  velocimeter  (LDV)  measurements 
were  taken,  providing  mean  velocities  and  Reynolds  stresses  at  various  a^al 
locations  downstream  of  the  flameholder.  The  CFD  analysis  consisted  of  2-D  and  3- 
D  time-accurate  simulations,  with  various  turbulence  models  including  standard  k- 
e  model,  RNG  k-e  model,  and  full  Reynolds  stress  model.  The  analysis  showed  that 
the  3-D  effects  from  the  presence  of  the  walls  were  significant.  3-D  time  accu^te 
calculations  with  the  RNG  k-e  turbulence  model  were  necessary  to  accurately  predict 
the  isothermal  flowfield.  Good  overall  agreement  was  shown  between  CFD  analysis 
and  measurements. 

A  water-cooled  combustion  test  section  was  designed  and  fabricated  to  facilitate 
collection  of  CARS  and  LDV  data  at  Wright  Laboratory  for  combusting  flameholder 
flows.  After  several  delays,  the  tests  were  canceled,  due  to  the  lack  of  availability  of 
the  test  facilities.  In  their  place,  further  CFD  studies  were  completed,  modeling  the 
reacting  flow  past  a  triangular  flameholder  in  a  duct  (experiment  by  Sjunnesson,  et 
al.).  The  results  showed  that  the  acoustic  interactions  between  the  combusting  flow 
and  the  duct  were  of  paramount  importance  at  low  (near  blowout)  equivalence 
ratios,  and  that  small  changes  in  the  overall  equivalence  ratio  dramaticdly  affected 
the  vortex  shedding  in  the  wake  of  the  flameholder.  The  CFD  analysis  showed  a 
longitudinal  combustion  instability  at  a  lower  equivalence  ratio  than  was  seen 
experimentally.  More  work  is  needed  to  accurately  predict  reacting  flow  in  ducts. 
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Overall,  the  program  was  successful  in  developing  a  useful  database  for  isothermal 
bluff  body  flows,  suitable  for  CFD  validation.  The  use  of  time-accurate,  3-D  CFD 
analysis  was  shown  to  be  necessary  to  accurately  capture  the  unsteady  flowfield 
downstream  of  the  flameholder.  In  addition,  the  program  showed  that  rambumer 
performance  gains  may  still  be  possible  using  vortex  generators  with  the 
flameholder,  but  that  smdler  gains  must  be  accepted  to  avoid  large  pressure  losses. 
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1.0  INTRODUCTION 


1.1  Background 

The  advent  of  high  speed  cruise  vehicles  has  brought  about  the  need  to  develop 
efficient  propulsion  systems  for  the  Mach  0-6  flight  regime.  The  turboramjet  is  one 
such  propulsion  system.  In  its  simplest  form,  the  turboramjet  is  a  gas  turbine  with  a 
bypass  ramjet  built  into  the  system.  From  take-off  conditions  through  Mach  4  flight, 
the  turbine  engine  powers  the  vehicle.  Above  Mach  3,  the  ramjet  takes  over  and 
powers  the  vehicle  up  to  Mach  6  or  higher. 

The  development  of  the  ramburner  has  several  challenges.  One  of  the 
technological  challenges  is  the  development  of  fuel  injection /flameholding 
concepts  capable  of  providing  good  fuel-air  mixing,  stable  combustion 
characteristics,  and  low  pressure  loss.  Conventional  ramburners  have  used  dump 
combustors,  where  the  fuel  is  injected  from  the  wall.  But  dump  combustors  often 
staffer  from  poor  fuel-air  mixing  due  to  a  lack  of  adequate  fuel  penetration,  and  poor 
performance  due  to  an  insufficient  number  of  flameholders.  Instream  fuel 
injection  and  flameholding  is  thus  preferred.  Gas  turbine  augmentors  employ 
instream  devices.  However,  the  conventional  spraybars  and  v-gutter  flameholders 
used  in  augmentors  are  not  viable  for  high  speed  flight  due  to  constraints  placed  on 
the  system  by  the  high  temperatures  associated  with  Mach  3+  flight.  At  high  Mach 
numbers,  the  air  entering  the  combustor  is  hot  enough  to  autoignite  the  fuel  at  the 
injector,  possibly  causing  structural  damage  to  the  downstream  flameholder. 
Therefore,  the  instream  fuel  injector  and  flameholder  must  be  close-coupled,  and 
actively  cooled. 

Due  to  these  requirements,  the  Integral  Fuel  injector /Flameholder  (IFF)  was  selected 
as  the  baseline  scheme  by  United  Technologies  Research  Center  (UTRC)  in  a  recent 
turboramjet  technology  program  for  the  Air  Force.^-^  A  schematic  of  the  IFF  is 
shown  in  Figure  1.  The  IFF  concept  was  originally  conceived  as  an  advanced 
augmenter  scheme  at  United  Technologies  Pratt  and  Whitney  for  the  PW1120 
engine,  and  UTRC  adapted  the  concept  for  ramburner  applications.  Since  the 
coupled  injection /flameholder  unit  is  cooled  by  the  incoming  fuel,  the  IFF  avoids 
structural  damage  caused  by  the  autoignition  of  the  fuel. 
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Figure  1.  Schematic  of  IFF  Concept 


While  the  IFF  is  a  good  first  design,  the  combustion  efficiency  of  the  IFF  is  generally 
mixing-limited.  Therefore,  concepts  with  enhanced  fuel-air  mixing  are  needed  to 
reduce  the  ramburner  length  and  increase  combustion  efficiency.  Fuel-air  mixing  is 
strongly  dependent  upon  the  fuel  penetration  into  the  airstream.  Low  fuel 
penetration  can  provide  good  flameholder  stability,  but  at  the  expense  of  fuel-air 
mixing.  Likewise,  high  fuel  penetration  can  result  in  improved  mixing,  but  creates 
lean  fuel  mixtures  directly  behind  the  flameholder,  leading  to  combustion 
instabilities. 

1.2  Flameholder  Theory 

Bluff  body  flameholders  have  been  studied  for  many  years,  due  to  their  practical 
application  of  stabilizing  flames  in  flowing  combustible  mixtures.  Flame 
stabilization  is  necessary  since  the  flowing  gas  stream  can  exceed  the  normal  flame 
speeds  of  the  fuel-air  mixtures  by  orders  of  magnitude.  The  purpose  of  the 
flameholder  is  to  produce  a  recirculation  zone  in  which  the  combusting  gases  can 
reside  and  provide  continuous  ignition  to  the  flowing  mixture. 

The  detailed  mechanisms  of  flameholding  have  been  studied  extensively,  by 
examining  the  flow  around  the  bluff  body  itself.  In  isothermal  conditions,  the  flow 
in  the  wake  of  a  bluff  body  consists  of  alternating  vortices  shed  from  the  top  and 
bottom  of  the  bluff  body.  The  vortices  are  shed  at  a  nominal  Strouhal  number 
(/H/U)  of  0.2,  producing  a  symmetrical  double-vortex  flow  pattern  in  the  mean  flow 
of  the  wake.  However,  in  combusting  flow,  the  heat  released  from  the  combustion 
process  reduces  the  instantaneous  vortex  formation,  and  produces  a  steady  pair  of 
counter-rotating  vortices  in  the  wake.  As  the  equivalence  ratio  decreases  in  the 


flame  zone,  the  heat  release  also  decreases.  As  the  equivalence  ratio  in  the 
recirculation  zone  approaches  lean  flammability  limits,  the  alternating  vortex 
structure  reappears,  leading  to  flame  instabilities  that  eventually  extinguish  the 
flame. 

Models  of  flame  stability  have  generally  followed  two  lines  of  reasoning.  The  first 
follows  the  idea  that  the  recirculation  zone  in  the  wake  of  the  bluff  body  acts  as  a 
well  stirred  reactor,  with  additional  combustion  products  supplied  by  the  flow 
outside  of  the  recirculation  zone.  When  the  rate  of  supply  of  the  fuel-air  mixture 
exceeds  the  rate  at  which  it  can  be  burned,  the  flame  gets  extinguished.  The  second 
line  of  reasoning  holds  that  the  combustion  takes  place  in  the  shear  layer  between 
the  hot  recirculation  zone,  and  the  fresh  fuel-air  mixture  of  the  outer  flow. 
Combustion  is  initiated  as  the  two  streams  mix  in  the  shear  layer.  If  the  gas  is  not 
ignited  by  the  time  the  mixture  within  the  shear  layer  has  reached  the  end  of  the 
recirculation  zone,  combustion  instabilities  result.  In  other  words,  the  residence 
time  of  the  entrained  mixture  in  the  shear  layer  should  exceed  the  reaction  time  for 
stable  combustion.io 

Research  into  flameholder  stability  indicates  that  the  aerodynamic  blockage  plays  a 
major  role.  The  larger  the  flameholder,  the  larger  the  recirculation  zone,  and  the 
better  the  flameholder  stability.  However,  large  blockages  also  create  higher 
pressure  loss,  reducing  the  thrust  performance  of  the  engine.  Therefore, 
flameholder  design  becomes  a  tradeoff  between  combustion  stability  and  total 
pressure  loss. 

Typically,  simple  two-dimensional  shapes  have  been  used  for  flameholders,  such  as 
flat  plates,  v-gutters,  and  circular  and  rectangular  cylinders.  In  some  case, 
axisymmetric  shapes  have  been  used,  such  as  discs  and  cones.  There  have  been  a 
limited  number  of  studies  on  irregular  shaped  flameholders.  Williams,  et.  al.^i, 
have  studied  elliptically  shaped,  gutter-type  flameholders  (Figure  2),  with  spanwise 
slots  cut  near  the  trailing  edge  of  the  flameholder  to  facilitate  gas  sampling  in  the 
wake.  They  noticed  that  the  presence  of  the  slots  increased  the  entrainment  of  fresh 
mixture,  and  increased  the  turbulence  in  the  recirculation  zone.  Nicholson  and 
Fieldi2  attempted  to  modify  the  wake  characteristics  of  a  flat  plate  flameholder  by 
introducing  a  splitter  plate  at  the  centerline  of  the  flameholder  base.  They  reported 
only  minor  differences  between  the  performance  of  the  original  flat  plate  and  that  of 
the  modified  plate.  Other  researchers  examining  similar  splitter  plate  geometries  1 3 
have  noticed  that  the  splitter  plate  reduces  the  turbulence  in  isothermal  conditions 
by  stabilizing  the  vortex  shedding  in  the  wake.  Such  differences  would  not  be 
expected  imder  combusting  conditions.  Also  in  an  effort  to  modify  the  wake  flow  of 
a  flameholder,  Faili4  attached  streamwise  tabs  to  the  edges  of  a  disc  flameholder. 
The  tabs  had  very  little  effect  on  the  aerodynamic  blockage  of  the  disc,  and  had  little 
effect  on  the  flameholder  stability.  Stwalley,  et.  al.15  have  examined  the  use  of  v- 
gutter  flameholders  with  slots  of  various  widths  cut  into  the  trailing  edge.  The 
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purpose  was  to  examine  the  flameholding  characteristics  of  the  irregular  external 
surfaces  of  a  damaged  aircraft.  They  found  that  removing  material  from  the  trailing 
edge  of  the  flameholder  in  a  random  fashion,  ^ways  resulted  in  decreased  flame 
stability.  As  more  material  was  removed,  and  the  shape  approached  the  original 
shape,  albeit  with  a  smaller  characteristic  dimension,  the  flame  stability  improved. 
The  pressure  losses  were  reduced  as  material  was  removed,  but  this  was  a  simple 
function  of  reduced  blockage.  Rizk  and  Lefebvre^^  have  also  noticed  in  their  studies 
that  single  sided  flameholders  (such  as  a  plate  mormted  on  a  wall)  tend  to  provide 
better  combustion  stability,  but  more  pressure  loss,  that  a  double  sided  instream 
flameholder.  This  is  due  to  the  stable  single  vortex  produced  behind  a  single  sided 
flameholder. 


Figure  2.  Elliptical  -  Gutter  Flameholder  Tested  by  Williams  et  al. 

In  general,  the  research  shows  that  the  stability  of  a  flameholder  is  not  significantly 
affected  unless  the  aerodynamics  of  the  wake  are  altered.  Alteration  by  irregular 
geometry  often  leads  to  increased  turbulence  in  the  wake  and  reduced  stability 
characteristics. 
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1.3  Improved  Concepts 


In  the  current  study,  new  flameholder  concepts  have  been  introduced  that  build  on 
the  original  IFF  design.  The  new  concepts  all  use  the  enhancement  of  large  scale 
axial  mixing  to  improve  the  fuel-air  mixing  of  the  IFF.  The  idea  comes  from  the  use 
of  forced  mixers  developed  for  turbofan  applications  as  shown  in  Figure  3.  The 
mixers  provide  large  scale  axial  vortices  to  enhance  the  mixing  between  the  two 
fluids  on  each  side  of  the  plate.  The  mixers  also  increase  the  total  area  of  the  shear 
layer  between  the  two  fluids  by  convoluting  the  wake.  When  applied  to  the  IFF,  the 
two  fluids  above  and  below  the  mixer  are  the  same,  but  there  are  regions  of  high 
fuel  concentration  that  need  to  be  mixed  with  the  rest  of  the  air. 

Lobed  mixers  have  received  the  attention  of  the  research  community  in  recent 
years.  McCormick  and  Bennett^^  have  shown  experimentally,  that  small  scale 
horseshoe  vortices  are  generated  in  the  trough  of  the  mixer  plate,  further  enhancing 
mixing.  They  also  showed  that  the  thrust  of  the  system  was  actually  improved  in 
addition  to  the  increased  mixing.  Abdolfadl  and  Sehrais  performed  experimental 
tests  on  lobed  mixers,  varying  the  length  and  scallop  of  the  mixer  trailing  edge. 
They  found  that  the  proper  scalloping  can  increase  mixing.  Peschke  and  McVeyt9 
examined  the  use  of  lobed  mixers  positioned  above  a  back-step  flameholder  to 
enhance  mixing  and  combustion.  They  found  that  the  lobed  mixers  increased 
mixing  through  the  combustion  region,  increasing  the  flame  front  angle  by  a  factor 
of  four.  No  pressure  measurements  were  reported.  A  numerical  study  was  initiated 
by  Malecki  and  Lord20,  who  found  that  Navier-Stokes  methods  were  able  to 
properly  capture  the  physics  behind  the  lobed  mixer,  by  comparing  to  experimental 
data.  Finally,  Elliot,  et  al.2i,  conducted  a  joint  numerical-experimental  test,  and 
showed  that  lower  lobe-height /lobe- wavelength  ratio  produced  stronger  vortices. 
Although  lobed  mixers  have  proven  capabilities  in  providing  increased  mixing, 
their  application  to  ramjet  combustors  is  not  clear,  due  to  the  tradeoff  between 
mixing  and  stability. 

The  Phase  I  study  of  this  program  evaluated  two  possible  advanced  concepts  using 
the  lobed  mixer  approach.  The  Convoluted  Lip  IFF  is  shown  in  Figure  4,  along  with 
the  Split-Notched  IFF.  The  results  in  Phase  I  showed  the  potential  for  increased 
mixing,  but  did  not  directly  examine  the  issues  of  total  pressure  loss  and  flame 
stability. 
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Lobed  Mixer 

Figure  3.  Axial  Vortidty  Mechanism  of  Increased  Fuel-Air  Mixing 
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SNIFF 


Figtire  4.  Concepts  Numerically  Tested  in  Phase  I 
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2.0  TECHNICAL  OBJECTIVES  AND  APPROACH 


The  overall  goal  of  this  project  was  to  develop  a  fuel  injector/flameholder  design 
that  would  provide  increased  fuel-air  mixing  with  acceptable  increases  in  the  cold 
pressure  loss.  In  addition,  the  program  sought  to  develop  an  experimental  data  base 
for  the  validation  of  CFD  software.  The  baseline  design  chosen  for  this  study  was 
the  IFF  configuration  of  United  Technologies,  shown  in  Figure  1. 

A  parallel  CFD/ experimental  approach  was  defined  for  the  Phase  n  program  with 
the  following  objectives: 

1)  perform  3-D  numerical  studies  of  the  baseline  design  and  advanced 
designs.  Then  analyze  the  advanced  designs  and  select  the  best 
candidate(s)  for  isothermal  testing; 

2)  design,  fabricate,  and  experimentally  test  the  baseline  design,  and  two 
advanced  designs,  at  Wright  Laboratory;  compare  the  data  to  numerical 
results  for  possible  numerical  modeling  improvements; 

3)  perform  additional  3-D  numerical  studies  based  on  experimental 
results,  to  obtain  final  candidate  for  combustion  tests; 

4)  design  and  fabricate  combustion  test  section  to  facilitate  combustion 
tests;  and 

5)  design,  fabricate,  and  test  baseline  model  and  final  advanced  design 
under  reacting  conditions.  Tests  include  LDV  and  CARS  data  under 
non-reacting  and  reacting  conditions. 

The  first  four  objectives  were  successfully  accomplished  during  this  Phase  n  effort. 
The  fifth  task  was  not  performed  due  to  scheduling  difficulties  at  Wright  Laboratory. 
In  its  place,  reacting  CFD  analysis  was  performed  and  compared  to  existing 
benchmark  data. 
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3.0  ISOTHERMAL  BASELINE  EXPERIMENTAL  TESTS 
3.1  Overview 


The  purpose  of  the  isothermal  experimental  tests  was  to  establish  a  baseline  set  of 
data  for  the  validation  of  the  CFD  effort,  and  to  test  a  selected  set  of  advanced 
designs  for  possible  further  study  in  the  combustion  tests.  The  initial  test  matrix 
called  for  the  extensive  use  of  simultaneous  3-component  LDV,  along  with  Mie 
scattering  test  to  determine  species  concentration  levels  during  injection  tests.  This 
test  program  was  very  ambitious  and  would  have  provided  a  wealth  of  velocity, 
turbulence,  mbdng,  and  atomization  data  for  both  liquid  and  gaseous  fuel  injection. 
However,  many  difficulties  were  encoimtered,  and  the  test  matrix  had  to  be  severely 
curtailed.  Only  the  baseline  IFF,  without  injection,  was  thoroughly  tested  using  the 
LDV.  No  injection  tests  were  completed  and  no  Mie  scattering  results  were 
obtained. 

The  data  taken  for  the  baseline  tests  are  presented  in  three  main  groups:  1)  the 
single  component  LDV  data,  2)  the  three-component  near-field  LDV  data,  and  3)  the 
three  component  far-field  LDV  data.  The  single-component  data  were  taken  at  the 
13  far  field  axial  stations,  defined  as  x/H  =  -5,  -3, 1, 1.5,  2,  3,  4,  6,  8, 10, 14, 18,  and  22. 
Single  component  data  were  taken  to  allow  axial  velocity  measurements  from  the 
top  of  the  tunnel  to  the  bottom  of  the  ttmnel,  providing  a  check  on  the  conservation 
of  mass  from  station  to  station,  and  lending  confidence  to  the  three-component 
measurements.  The  three-component  measurements  were  taken  in  the  far  field, 
and  in  the  near  field,  defined  as  axial  stations  x/H  =  0.5,  0.6,  0.7,  0.8,  0.9,  1.0, 1.1,  and 
1.2. 

In  addition  to  the  baseline  tests,  one  advanced  model  (ADV2)  was  tested  using  the 
LDV  system.  Vertical  profiles  were  taken  at  x/H  =  2.5,  and  z/H  =  2.7,  2.0, 1.0, 0.0,  -1.0, 
and  -2.0.  Spanwise  profiles  were  taken  at  x/H  =  4.0,  and  y/H  =  1.5,  0.0,  and  -1.5. 
Several  other  model  variations  were  briefly  tested,  in  an  effort  to  determine  the 
peak  crossflow  velocity  in  the  near  wake  of  the  flameholder.  The  conclusions  from 
these  tests  are  presented,  along  with  an  analysis  of  the  measured  pressure  losses. 

3J2  Test  Facility 

The  test  facility  was  located  in  Test  Cell  18  at  the  Aero  Propulsion  and  Power 
Directorate  of  Wright  Laboratory.  The  test  stand  was  specifically  designed  and 
constructed  to  accommodate  flow  visualization  and  laser-based  diagnostics.  The  76 
mm  X  130  mm  (3.0  in  x  5.0  in)  test  section  originally  had  full-view  fused  quartz 
windows  on  both  sides,  with  a  partial  view  window  on  the  top.  For  the  current 
study,  the  bottom  of  the  tunnel  was  modified  to  also  include  a  partial  view  window. 
The  tunnel  air  was  supplied  by  the  facility  compressors,  allowing  total  pressure 
beyond  500  psia  at  ambient  temperature.  With  choke  blocks  located  downstream  of 
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the  test  section,  the  free  stream  Mach  number  was  fully  variable  between  0  and  1. 
Two  sets  of  honeycomb  and  screens  were  used  in  a  small  settling  chamber  upstream 
of  the  test  section  to  reduce  the  freestream  turbulence,  and  provide  uniform  flow  to 
the  test  section.  Facility  instrumentation  included  an  array  of  absolute  pressure 
transducers  and  thermocouples  monitored  by  a  micro-computer  based  system. 

The  original  test  conditions  were  designed  to  follow  closely  the  earlier  experimental 
tests  of  Hautman,  et  al.9,  but  at  ambient  temperature  (M  =  0.2,  p  =  275  BCPa,  T  =  300 
K).  However,  initial  testing  indicated  the  presence  of  a  low  frequency  pressure 
oscillation  in  the  test  section.  Much  time  and  effort  was  spent  attempting  to 
alleviate  the  problem.  A  new,  larger  settling  chamber  was  designed  and  fabricated  to 
damp  the  oscillations.  The  chamber  included  a  high  pressure  loss  core  buster, 
intended  to  act  as  an  acoustic  barrier  to  small  changes  in  the  supply  pressure. 
Unfortunately,  this  did  not  solve  the  problem. 

The  next  tunnel  modification  involved  powering  the  tunnel  from  an  auxiliary  air 
supply.  Several  options  were  explored,  including  the  use  of  a  Roots  blower  to 
provide  high  pressure  air.  However,  the  tunnel  was  finally  connected  to  a  7.5  Hp 
centrifugal  blower.  This  provided  air  pressures  only  slightly  above  ambient 
pressure  at  a  peak  Mach  number  of  0.08.  It  was  determined  that  the  lower  Mach 
number  would  be  sufficient  to  provide  baseline  CFD  data,  since  the  original  flow 
conditions  were  well  in  the  incompressible  regime.  The  problem  with  the  pressure 
oscillations  were  still  present,  however,  even  with  the  blower  supplying  the  air. 
Additional  testing  conduded  that  using  either  the  original  small  settling  chamber  or 
the  larger  one  did  not  drastically  affect  the  osdllations. 

Therefore,  attention  turned  to  the  components  of  the  tunnel  downstream  of  the  test 
section.  Since  the  choke  blocks  were  no  longer  being  utilized,  the  test  section  exit 
was  no  longer  acoustically  isolated  from  the  downstream  sections.  The  trials 
included:  1)  removing  the  exit  diffuser  and  allowing  the  test  section  to  exhaust 
directly  into  the  room,  2)  allowing  the  test  section  to  exhaust  into  a  dump  diffuser 
with  no  suction  applied,  and  3)  attaching  a  long  7  inch  diameter  flex  hose  to  the  test 
section  exit  and  exhausting  the  flow  through  a  powered  roof  vent.  The  last 
configuration  provided  the  most  stable  flow  conditions  although  they  were  not 
ideal.  The  final  configuration  chosen  for  the  tests  were  the  7.5  Hp  blower  with 
honeycomb  at  the  blower  exit,  ducted  into  the  small  settling  chamber,  and  into  the 
test  section  through  another  layer  of  honeycomb. 

For  the  final  tunnel  configuration,  the  free  stream  dynamic  pressure  was  measured 
8  inches  upstream  of  the  test  section  with  a  pitot  static  probe  connected  to  both  sides 
of  a  micromanometer.  The  total  temperature  was  measured  at  the  same  axial 
location  with  a  thermocouple  inserted  in  the  flow,  and  the  static  pressure  was 
measured  with  an  absolute  pressure  transducer  connected  to  the  side  wall  of  the 
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tunnel.  All  reference  conditions  were  recorded  by  hand  during  the  acquisition  of 
each  data  point  in  the  test  matrix. 

- Model  Definition  - - 

The  baseline  IFF  model  was  fabricated  out  of  aluminum  in  three  sections  (Figure  5). 
The  hemispherical  nose  had  a  radius  of  0.25  inches.  The  rectangular  center  section 
followed  with  a  length  of  0.625  inches,  and  housed  a  fuel  manifold  for  the  injection 
tests.  Finally  a  detachable  base  plate  measuring  0.125  inches  long  was  connected  to 
the  trailing  edge.  The  base  plate  was  detachable  to  allow  the  attachment  of  modified 
trailing  edges  to  the  IFF.  The  model  was  mounted  in  the  test  section  directly  to  the 
side  wall  windows  4  inches  downstream  of  the  test  section  inlet,  allowing  free 
optical  access  to  the  entire  geometry.  The  reference  point  for  all  data  is  defined  as 
the  vertical  and  spanwise  centerlines  at  the  IFF  trailing  edge.  The  total  length  of  the 
IFF  was  1.0  inches  with  a  base  height  H  =  0.5  inches,  giving  a  total  blockage  of  17%  in 
the  3  inch  high  test  section. 


not  to  scale 


fuel  manifold 


R  =  0.25"' 


pm  connection 


screw  connection 


0.5" 


Figure  5.  Experimental  Model  Design 
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3.4  Inlet  Conditions 


The  nieasured  inlet  conditions  changed  slightly  from  test  to  test,  since  the  entire  test 
spanned  a  period  of  ol  year.  Each  data  point  was  nondimensionaiized  with  the 
actual  freestream  conditions.  For  reference  purposes  the  following  inlet  conditions 
were  specified. 

Uref  =  24.5  m/sec 
Tref  =  295K 
Ren  =  18,000 
Mref  =  0.07 


3.5  LDV  Hardware 

Two  different  optical  systems  were  used  for  the  far-field  and  the  near-field  tests.  The 
differences  were  all  contained  within  the  transmitting  optics;  the  set  up  for  the 
receiving  optics  was  the  same  for  all  tests.  The  entire  system  was  mounted  on  a 
three-axis  motion  table  capable  of  placing  the  measurement  probe  volume  at  any 
point  in  the  test  section  with  a  resolution  of  ±  0.03  mm. 

3.5.1  Fat-Field  Transmitting  Optics 

The  far  field  LDV  system  started  with  a  1  watt  argon-ion  laser.  The  laser  beam  was 
steered  into  a  TSI  Color  Burst  beam  separator,  where  the  beams  were  separated  into 
their  primary  colors,  polarized,  split,  and  Bragg  shifted.  Fiber  optic  couplers 
transmitted  the  conditioned  beams  to  the  transmitting  heads  adjacent  to  the  test 
section.  The  blue  and  green  beams  were  directed  through  the  side  windows  from  a 
TSI  model  9277  two-component  fiber  optic  probe  with  a  focal  length  of  450  mm, 
such  that  the  green  beams  ran  parallel  to  the  tunnel  and  the  blue  beams  ran 
perpendicular  to  the  tunnel.  The  violet  beams  were  directed  through  the  top 
window  via  a  TSI  model  9274  single-component  fiber  optic  probe.  All  three 
components  were  Bragg  shifted  40  MHz  to  resolve  the  directional  bias. 

To  ensure  that  all  six  beams  from  the  three  different  components  were  focused  at 
the  same  point  in  space,  the  beams  were  steered  to  pass  through  a  35  |im  pinhole 
oriented  at  45  degrees  to  the  vertical.  The  alignment  was  then  further  optimized  by 
maximizing  the  3-component  data  rate  from  a  w’ater  mist  passing  through  the  probe 
volume. 

3.5.2  Neai-Field  Transmitting  Optics 

During  the  time  between  the  far-field  tests  and  the  near-field  tests,  the  two- 
component  fiber  optic  probe  was  damaged.  Although  the  unit  was  repaired,  it  was 
not  returned  in  time  for  the  near  field  tests.  Therefore,  the  near  field  tests  had  to 
use  an  alternative  optical  train  for  the  transmitting  optics.  The  alternative  method 
involved  putting  together  a  manual  optic  train  from  components  available 


14 


throughout  the  laboratory.  The  process  was  very  tedious  but  ultimately  successful. 
Figure  6  shows  a  schematic  of  the  system  used  for  the  near  field  tests.  The  beam 
from  the  argon-ion  laser  was  first  separated  into  the  primary  colors  with  a  prism. 
The  blue  and  green  beams  were  then  processed  separately  through  a  polarizer,  beam 
splitter,  and  Bragg  cell,  before  being  focused  into  the  test  section.  The  beam 
orientation  was  rotated  45  degrees  for  the  near  field  set  up  allowing  the  probe 
volume  closer  access  to  the  base  of  the  IFF.  The  violet  beam  from  this  laser  was  not 
used. 

A  second  1  watt  argon-ion  laser  was  used  to  generate  the  violet  beam  in  the  single 
color  mode.  This  provided  for  a  slight  increase  in  the  power  available  to  the  violet 
beam,  which  was  found  to  be  a  limiting  factor  in  the  far  field  tests.  The  violet  beam 
was  then  passed  through  the  color  burst  separator  to  take  advantage  of  the  bragg  cell 
internal  to  the  system.  The  violet  beams  were  then  connected  via  fiber  optic  cables 
to  the  single  component  fiber  optic  probe  and  utilized  as  in  the  far  field  tests.  Beam 
alignment  between  the  different  colors  was  accomplished  in  a  manner  similar  to  the 
far  field  tests,  and  all  components  were  Bragg  shifted  40  MHz. 


Figme  6.  Near-Field  Optical  Set-Up  for  3-D  LDV  System 
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3.5.3  LDV  Receiving  Optics 

All  LDV  channels  were  operated  in  a  forward  scattering  mode  to  maximize  the 
collected  signals.  The  blue  and  green  signals  were  collected  with  a  120  mm  lens  and 
passed  through  a  dichroic  filter.  The  filter  passed  the  green  light  into  a 
photomultiplier,  and  reflected  the  blue  light  into  another  photomultiplier.  The 
violet  light  was  also  collected  in  the  forward  scattering  mode  and  passed  to  a  third 
photomultiplier. 

The  output  from  each  photomultiplier  was  routed  to  TSI's  IFA  750  Digital  Burst 
Correlator  and  Analyzer.  The  IFA  750  was  able  to  collect  the  data  at  a  high  rate  of 
speed  while  performing  a  number  of  validation  checks  on  each  velocity  realization. 
T^e  signals  were  first  conditioned  by  passing  through  high  and  low  pass  filters,  set  at 
20  MHz  and  100  MHz  respectively,  for  each  channel.  One  of  the  validations 
included  the  use  of  a  user  controlled  coincidence  window,  defined  as  the  amount  of 
time  allowed  for  each  channel  to  produce  a  valid  signal  for  a  single  particle. 
Realizations  from  different  channels  outside  of  the  coincidence  window  are  deemed 
to  be  from  different  particles  and  discarded.  A  window  of  20  usee  was  selected  for  all 
the  current  LDV  tests. 

Another  check  was  a  count  of  the  number  of  fringes  crossed  by  a  given  particle.  The 
user  could  set  a  minimum  number  of  fringes  that  the  particle  must  cross,  or  the 
signal  was  invalid.  If  any  tests  failed  for  any  one  of  the  channels,  then  the  point  was 
invalidated  for  all  channels.  A  minimum  of  8  fringes  was  chosen  for  this  study. 
Valid  data  points  were  collected  in  a  data  buffer  within  the  IFA  750  and  downloaded 
to  the  controlling  computer  in  blocks  of  1024  points. 

3.5.4  Data  Acquisition  and  Postprocessing 

A  486-based  personal  computer  was  used  to  control  the  EFA  750  through  the  use  of 
TSI's  FIND  software  package.  The  software  allowed  the  user  to  control  every  aspect 
of  the  IFA  750,  and  to  view  real-time  histograms  of  the  incoming  signals.  After  the 
data  was  acquired,  the  software  provided  full  data  reduction  and  plotting 
capabilities. 

During  the  course  of  the  testing,  some  minor  shortcomings  were  encountered  with 
both  the  acquisition  and  reduction  portions  of  the  software.  During  data  acquisition, 
the  coincidence  window  was  not  always  set  properly,  even  though  the  software 
indicated  that  it  was  set.  Cycling  the  coincidence  window  setting  seemed  to  solve 
the  problem.  The  difficulties  with  the  reduction  software  were  more  difficult  to 
alleviate.  The  FIND  software  made  several  assumptions  that  did  not  hold  with  the 
LDV  configuration  used  in  the  current  study.  First,  the  software  assumed  that  the 
green  and  blue  components  were  perfectly  orthogonal.  Although  this  should  be 
true,  measurements  indicated  that  the  beams  were  nonorthogonal  by  approximately 
1  degree.  Second,  the  software  assumed  that  any  frequency  shifting  would  be  in  the 
negative  direction,  relative  to  the  local  coordinate  system.  While  this  is  generally 
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true  for  low  speed  axial  measurements,  it  is  arbitrary  for  flows  with  a  net  zero 
velocity,  such  as  the  cross-flow  velocity  components  in  a  duct.  In  the  current 
configuration,  the  blue  beam  was  shifted  in  the  positive  direction.  The  FIND 
software  would  not  accept  a  negative  frequency  shift,  but  would  accept  a  negative 
fringe  spacing.  Such  a  solution  is  mathematically  copect,  however  TSI 
recommended  against  using  a  negative  fringe  spacing,  indicating  that  the  software 
often  uses  hidden  negative  signs  on  certain  variables  to  flag  unrelated  operations. 

To  resolve  the  reduction  software  problems,  a  new  data  reduction  package  was 
written  to  run  independent  of  the  TSI  software.  The  new  software  was  developed 
on  a  486  based  PC  using  the  Lahey  FORTRAN  compiler,  and  read  the  raw  data 
directly  from  the  TSI  binary  packed  files.  The  software  included  full  control  over  all 
beam  angles,  and  frequency  shift  directions  and  magnitudes.  It  also  allowed  for  the 
inclusion  of  reference  conditions  to  facilitate  point-by-point  normalization  of  the 
data  The  new  software  was  tested  and  validated  against  the  FIND  software  and 
produced  identical  results  for  the  test  cases.  Additional  information  on  the  new 
software,  including  the  source  code,  is  presented  in  Appendix  C. 

During  the  calculation  of  the  flow  statistics,  additional  measures  were  taken  to 
remove  spurious  points.  The  velocity  measurements  were  first  transformed  to  the 
tunnel  coordinate  system  for  each  validated  realization.  The  Probability  Density 
Functions  (PDF)  were  then  assembled  for  each  velocity  component  at  the  given 
measurement  stadon.  The  statistics  of  interest  were  calculated  by  the  methods 
described  in  Appendix  B.  Finally,  all  velocity  points  were  compared  against  a  3 
sigma  criteria.  Any  points  falling  outside  of  3  sigma  from  the  mean  were  discarded. 
If  a  point  failed  for  any  one  of  the  three  velocity  components,  then  the  point  was 
removed  from  all  components.  All  statistics  were  then  recomputed  with  the 
spurious  points  removed.  For  the  current  study,  5120  validated  points  were 
collected  from  the  IFA  750  at  each  station.  Of  these,  typically  50  points  were 
eliminated  by  the  3  sigma  criterion. 

Another  source  of  error  for  LDV  measurements  is  velocity  bias.  Velocity  bias  arises 
from  the  fact  that  turbulent  flow  is  continuous  and  structured.  Therefore,  for  a 
uniformly  seeded  flow,  during  the  period  of  time  that  the  flow  is  traveling  faster 
than  the  true  mean,  the  LDV  will  measure  more  particles  than  when  the  flow  is 
traveling  slower  than  the  true  mean.  The  simple  arithmetic  mean  calculated  from 
the  collected  measurements  will  then  be  artificially  high.  The  higher  the  turbulence 
level,  the  greater  the  effect  of  the  bias.  Although  several  methods  have  been 
devised  to  correct  for  velocity  bias,  general  acceptance  of  the  post  processing 
corrections  is  lacking.  While  the  FIND  software  included  velocity  bias  corrections, 
the  new  reduction  software  did  not  include  such  corrections,  and  no  velocity  bias 
corrections  were  made  to  the  current  data  set. 
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3.5.5  LDV  Seed  Selection 

The  next  problem  involved  the  proper  selection  of  seed  material  to  use  with  the 
LDV  system.  Several  different  seed  materials  were  tested.  Initial  tests  used  the 
chemical  seeder  designed  by  Craig,  et.  al.22  The  chemical  seeder  produced  small 
titanium  dioxide  particles  (0.2  -  1  )im)  by  reacting  titanium  tetrachloride  with  the 
water  in  the  air  supply,  according  to  the  following  equation: 

Tick  +  2(H20)  ->  Ti02  +  4(HCl) 

The  titanium  dioxide  particles  produced  by  the  reaction  make  excellent  seed 
particles,  but  titanium  tetrachloride,  and  the  reaction  by-product,  hydrochloric  acid, 
are  extremely  corrosive  and  difficult  to  handle.  Due  to  the  hazards  of  using  such 
corrosive  seed  material,  alternative  seeding  methods  were  explored. 

Aluminum  oxide  powder  was  the  next  seed  material  used  for  the  tests.  This  finely 
ground  powder  supplied  seed  particles  from  0.5  -  5  microns  in  diameter,  provided 
the  material  was  kept  very  dry.  A  new  seeding  device  was  fabricated  using  a  hot 
plate  and  dry  nitrogen  purge  system  to  keep  the  seed  dry.  The  nitrogen  system  was 
then  used  to  create  a  fluidized  bed  to  inject  the  seed  into  the  tunnel.  While  the 
system  worked  well  at  times,  it  was  prone  to  clogging  and  would  often  surge, 
injecting  large  amounts  of  seed  into  the  tunnel.  Since  uniform  seed  rates  are 
desirable  during  data  sampling,  and  from  sample  to  sample,  the  surging  problem 
was  unacceptable.  The  aluminum  oxide  seed  also  had  a  tendency  to  build  up  on  the 
IFF  leading  edge  and  would  quickly  coat  the  windows  rendering  the  LDV  system 
useless. 

The  last  seeder  tested  was  of  the  liquid  glycerin  type,  and  was  fabricated  by  Rabe  and 
Sabroske.23  The  liquid  glycerin  seeder  produced  liquid  droplets  of  glycerin  that 
ranged  from  0.5  -  1.0  |im  in  diameter.  The  seed  rate  was  easily  controlled  and 
consistent.  In  general,  the  seed  material  did  not  collect  on  the  windows,  although 
the  separation  region  directly  behind  the  IFF  became  coated  fairly  quickly.  When 
taking  measurements  in  this  area,  the  windows  had  to  be  cleaned  after  every  3-4 
data  points. 

Peak  data  rates  with  all  the  seeders  were  between  15,000  and  20,000  per  second  in  the 
single  channel  mode.  During  simultaneous  3-component  data  collection,  the  data 
rate  dropped  to  5000  per  second  in  the  free  stream,  and  1500  per  second  in  the  highly 
turbulent  region  behind  the  IFF.  The  windows  were  removed  and  cleaned 
whenever  the  data  rate  fell  below  1000  per  second. 
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3.6  Experimental  Results  -  Baseline  IFF 

The  data  were  taken  in  three  main  groups,  designated  as  the  one-component  far 
field  data,  the  three  component  far  field  data,  and  the  three  component  near  field 
data.  The  far  field  data  were  obtained  with  each  color  aligned  directly  with  the 
velocity  component  being  measured.  Therefore,  no  coordinate  transformation  was 
required.  However,  such  an  arrangement  prevented  the  measurement  of  normal 
velocities  within  0.5  inches  of  a  wall,  due  to  partial  blockage  of  the  incoming  laser 
beams.  The  near  field  data  were  taken  with  the  blue  and  green  beams  rotated  45 
degrees  to  the  tunnel  axis,  to  allow  measurements  to  within  0.25  inches  of  the  IFF. 
The  data  then  had  to  be  transformed  to  the  tuimel  coordinate  system  for  processing. 

3.6.1  One-Component  Far  Field  Data 

The  first  data  collected  were  the  axial  component  at  the  far  field  stations.  Since  only 
the  axial  component  was  collected,  data  were  taken  in  profiles  from  the  floor  to  the 
ceiling  of  the  test  section.  The  velocity  and  turbulence  profiles  are  presented  in 
Figures  7  and  8,  and  show  several  basic  features  of  the  flow.  First,  the  incoming  flow 
profile  is  relatively  flat  and  symmetric,  indicating  uniform  inlet  conditions.  Second, 
the  freestream  turbulence  is  approximately  3%,  with  a  maximum  of  8%  in  the 
tunnel  wall  boundary  layers.  Finally,  the  presence  of  the  wake  is  clearly  seen  in  the 
velocity  deficit  and  increased  turbulence  levels  downstream  of  the  IFF.  Since  no 
negative  mean  velocities  were  indicated  in  the  wake,  the  recirculation  zone  ended 
within  IH  of  the  IFF. 

From  the  axial  velocity  profiles,  the  mass  flux  at  each  station  was  integrated  by  the 
trapezoidal  rule,  assuming  that  the  flow  was  2-D  with  constant  density.  The  mass 
flow  at  each  station  agreed  to  within  2%  of  the  first  station,  with  the  highest  errors 
-occurring  in  the  near  wake  region  (x/H  =  1,  1.5,  2).  Such  a  result  gives  a  indication 
that  the  measurements  were  taken  correctly,  and  that  velocity  bias  was  not  a 
problem.  However,  velocity  bias  will  only  be  significant  in  regions  of  high 
turbulence,  such  as  directly  behind  the  IFF.  But  the  region  in  the  near  wake  also 
makes  the  smallest  contribution  to  the  mass  flow,  due  to  the  velocity  deficit. 
Therefore,  velocity  bias  may  still  be  significant  even  though  the  mass  flow  errors  are 
small. 
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)(/H=10  x/H=14  x/H=18  x/H=22 


Several  samples  along  the  centerline  were  analyzed  using  Fourier  analyses  to 
determine  the  vortex  shedding  frequency  behind  the  IFF.  Since  the  LDV  data  was 
taken  by  random  samples  whenever  a  particle  passed  through  the  probe  volume,  a 
transformation  to  uniform  time  intervals  between  data  realizations  was  necessary. 
The  FIND  software  provided  a  mechanism  that  averaged  all  data  within  a  given 
time  interval,  dt.  The  entire  data  set  was  averaged  to  give  one  realization  for  every 
dt  within  the  data  set.  A  Fast  Fourier  Transform  (FFT)  routine  then  computed  the 
frequency  spectrum.  Figure  9  shows  the  frequency  spectrum  as  measured  at  x/H  =  -5 
and  2.  A  distinct  frequency  peak  is  noted  at  454  Hz  (St  =  0.243)  for  the  downstream 
location,  while  the  upstream  location  shows  no  dominant  frequency.  The 
dominant  frequency  was  found  to  continue  at  all  downstream  stations  at 
diminishing  magnitudes  imtil  it  disappeared  at  the  last  station,  x/H  =  22. 

3.6.2  Three-Component  Far  Field  Data 

After  the  single-component  data  were  collected,  all  three  components  were 
collected  at  the  far  field  stations.  The  general  agreement  for  the  axial  component 
between  the  single  component  and  three-component  measurements  was  very  good. 
A  comparison  is  shown  in  Figure  10  at  the  near  wake  station  of  x/H  =  1.0.  This 
station  is  the  most  demanding  station  of  the  far  field  data,  due  to  the  high 
turbulence,  and  there  is  some  discrepancy  between  the  two  data  sets.  While  several 
possible  errors  were  identified,  the  greatest  error  was  due  to  the  replacement  of  the 
model  in  the  test  section  between  the  two  tests.  Every  effort  was  made  to  maintain 
the  model's  attitude  and  position  within  the  test  section,  but  small  errors  were 
possible. 

In  addition  to  the  conclusions  drawn  from  the  single-component  data,  the  three- 
component  data  (presented  in  Figures  11-16)  yielded  the  following  information. 
The  vertical  velocity  profiles  indicated  that  the  flow  is  turning  back  towards  the 
centerline  in  the  wake  as  the  flow  fills  the  low  pressure  region  left  behind  the  IFF. 
The  maximum  velocity  gradients  in  the  vertical  direction  were  located  at  the 
centerline,  and  also  marked  the  position  of  the  highest  turbulence  intensity  (75%  at 
x/H  =  1).  Likewise,  the  axial  component's  maximum  velocity  gradient  was  located 
in  the  shear  layer  at  y/H  =  ±0.2,  where  the  turbulence  was  also  at  its  peak.  The 
regions  of  maximum  turbulence  in  the  axial  and  transverse  components  did  not 
spatially  correspond,  indicating  that  the  turbulence  was  not  isotropic. 
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Figure  9.  Frequency  Spectra  Upstream  and  Downstream  of  IFF  Shows  Vortex  Shedding  Frequency 
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Normalized  Mean  Transverse  Velocity  Profiles  at  Far  Field  Stations;  3-D  LDV 
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The  reason  the  turbiiience  was  not  isotropic  is  due  to  the  structure  inherent  in  the 
turbulence.  The  turbulence  in  the  IFF  wake  is  dominated  by  large  scale  vortices 
shed  alternately  from  the  top  and  bottom  comers  of  the  IFF  trailing  edge.  As  each 
vortex  is  shed,  they  form  a  train  of  vortices  of  alternating  directions  moving 
downstream  through  the  wake.  The  turbulence  is  now  simply  a  measure  of  the 
velocity  variations  due  to  the  vortices  moving  through  a  sampling  point.  If  for 
example,  the  sampling  point  is  on  the  centerline,  there  will  be  very  little 
contribution  to  the  axial  component  from  the  vortices,  since  the  vortex  cores  are 
traveling  near  the  centerline.  However,  there  is  a  large  contribution  to  the  vertical 
velocity.  The  vertical  velocity  varies  between  ±Vmax  as  each  vortex  passes  the 
sampling  point,  leading  to  high  turbulence  at  the  centerline.  This  also  leads  to  a 
bimodal  PDF  such  as  the  one  measured  at  x/H  =  1,  y/H  =  0  (Figure  17).  Likewise,  the 
highest  variation  in  the  axial  component  occurs  at  the  edges  of  the  mean  shear 
layer,  where  very  little  vertical  velocity  variation  exists. 

The  spanwise  velocity  component  (W)  was  negligible  throughout  most  of  the  test 
measurement  domain,  although  a  small  mean  flow  was  measured  at  x/H  =  1.  Two- 
dimensional  geometries  have  been  observed  to  produce  three-dimensional  flows 
due  to  the  presence  of  the  tunnel  boundary  layers.  However,  the  spanwise 
centerline  velocity  should  have  no  mean  spanwise  velocity  component.  The  fact 
that  there  was  spanwise  flow  measured  at  x/H  =  1  indicates  one  of  two  possibilities: 
1)  the  flow  and/or  the  geometry  was  not  exactly  symmetric,  or  2)  the  z  =  0  station 
was  not  exactly  centered  with  respect  to  the  flow.  TTie  location  of  the  reference  point 
was  accomplished  by  focusing  the  probe  volume  on  each  side  wall  and  splitting  the 
difference.  The  source  of  the  spanwise  flow  was  further  investigated  in  the  near 
field. 

3.6.3  Three-Component  Near  Field  Data 

The  near  field  three-component  data  are  presented  in  Figures  18-26.  The  mean 
velocity  profiles  indicated  that  the  measurement  stations  are  within  the 
recirculation  zone,  and  that  the  recirculation  zone  extended  to  x/H  =  0.9.  The  length 
of  the  recirculation  zone  is  much  smaller  than  that  of  a  backward  facing  step,  due  to 
the  instantaneous  structure  of  the  vortex  shedding  process  that  is  not  present  with  a 
backward  fadng  step. 

The  profiles  of  the  Reynolds  normal  stresses  show  that  the  peak  intensities  occurred 
at  locations  where  the  mean-velocity  gradients  were  the  highest.  The  normal 
stresses  are  also  clearly  not  isotropic,  as  seen  in  a  comparison  of  the  stress  ratios  in 
Figures  27-29.  The  vertical  normal  stress  (vv)  was  found  to  be  the  dominate  normal 
stress  term.  Reynolds  shear  stresses  are  also  shown.  The  measurements  indicated 
that  the  uv  shear  stress  is  the  most  significant  of  the  three,  and  is  five  times  higher 
than  the  other  two  shear  stresses. 
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Normalized  Mean  Axial  Velocity  Profiles  at  Near  Field  Stations;  3-D  LDV 
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Figure  19.  Normalized  Mean  Transverse  Velocity  Profiles  at  Near  Field  Stations;  3-D  LDV 
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Normalized  Mean  Spanwise  Velocity  Profiles  at  Near  Field  Stations;  3-D  LDV 
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Figure  22.  Normalized  vv  Turbulent  Stress  Distributions  at  Near  Field  Stations 


38 


39 


Normalized  uv  Turbulent  Stress  Distributions  at  Near  Field  Stations 


Figure  25.  Normalized  uvv  Turbulent  Stress  Distributions  at  Near  Field  Stations 
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Figure  27.  Axial  Mean  Normal  Stresses  -  Isotropy  Check  (Isotropic  uu/k  =  2/3) 


x/H=0.5  x/H=0.6  x/H=0.7  x/H=0.8 


Transverse  Mean  Normal  Stresses  -  Isotropy  Check  (Isotropic  vv/k  =  2/3) 


Also  present  in  the  near  field  data  is  the  spanwise  flow  in  the  wake  of  the  IFF.  It  was 
hypothesized  that  the  presence  of  the  tunnel  boundary  layers  created  a  pressure 
force  towards  the  center  of  the  tunnel  as  the  flow  passed  over  the  IFF.  The  flow 
reacts  to  the  pressure  force  by  forming  four  coxinter-rotating  vortices  on  the  back  face 
of  the  IFF,  as  viewed  from  downstream.  To  test  the  hypothesis,  several  one- 
component  (W)  measurements  were  made  in  an  attempt  to  capture  the  presence  of 
the  vortices  (Figure  30).  The  flow  can  clearly  be  seen  moving  towards  the  centerline 
at  the  very  near  stations,  and  away  from  the  centerline  downstream  of  these 
stations,  indicating  the  presence  of  at  least  two  counter-rotating  vortices,  one  on 
each  side  of  the  spanwise  centerline.  In  addition,  the  vertical  profiles  show  that 
there  are  two  counter-rotating  vortices  in  the  vertical  direction,  with  one  on  each 
side  of  the  vertical  centerline.  Therefore,  the  original  hypothesis  is  supported  by  the 
data. 
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4.0  ISOTHERMAL  NUMERICAL  INVESTIGATION  OF  BASELINE  IFF 
CONFIGURATION 

The  isothermal  test  matrix  for  the  baseline  IFF  can  be  divided  into  seven 
incompressible  cases.  The  first  two  cases  were  completed  assuming  that  the 
flowfield  was  2-D  and  symmetric;  only  the  top  half  of  the  geometry  was  modeled 
and  the  solution  was  converged  to  steady  state.  The  next  three  cases  modeled  the 
entire  geometry  in  2-D,  and  utilized  time-accurate  techniques  to  capture  the 
unsteady  flowfield  behind  the  IFF.  The  last  2  cases  were  3-D  solutions.  Table  1 
summarizes  the  isothermal  cases  analyzed. 

Table  1.  Isothermal  Numerical  Test  Cases 


Case 

Dimension 

Grid 

Turbulence  Model 

Time  Integration 

I 

2-D 

half-height 

k-e 

steady  state 

n 

2-D 

half-height 

full  Reynolds  stress 

steady  state 

m 

2-D 

full 

laminar 

time  accurate 

IV 

2-D 

full 

k-e 

time  accurate 

V 

2-D 

full 

RNG  k-e 

time  accurate 

VI 

3-D 

full 

RNG  k-e 

time  accurate 

VII 

3-D 

full 

RNG  k-e 

time  accurate 

All  solutions  used  2nd  order  accurate  central  differencing,  stabilized  with  upwind 
biasing.  The  time-accurate  solutions  were  integrated  through  time  using  the  Crank- 
Nicholson  technique. 

4.1  CFD  Code 

The  numerical  analyses  for  the  current  study  was  completed  using  a  commercially 
available  CFD  code,  CFD-ACE.24  Developed  by  CFD  Research  Corporation,  CFD-ACE 
includes  the  following  features: 

a.  solution  of  two-  and  three-dimensional  Navier-Stokes  equations  for 
incompressible  and  compressible  flows; 

b.  cartesian  and  nonorthogonal  curvilinear  coordinates; 

c.  single  and  multiblock  grid  topology; 

d.  steady-state  and  time  accurate  solution  algorithms; 

e.  advanced  turbulence  modeling,  including  standard  k-e,  RNG  k-e,  and 
full  Reynolds  stress  models; 

f.  instantaneous,  one-step,  two-step  and  four-step  combustion  models 
with  equilibrium  products  option,  and  prescribed  and  Monte  Carlo  PDF 
combustion  models; 
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g.  first  and  second  order  upwind,  central  with  damping,  and  Osher- 
Chakravarthy  differencing  schemes; 

h.  co-located,  fully  implicit  and  strongly  conservative  finite  volume 
formulation; 

i.  pressure-based  solution  algorithms  including  SIMPLE  and  a  variant  of 
SIMPLEC;  and 

j.  modified  form  of  Stone's  strongly  implicit  solver  and  conjugate 
gradient  solver. 


CFD-ACE  has  undergone  a  considerable  amount  of  systematic  quantitative 
validation  for  both  incompressible  and  compressible  flows.  Over  50  validation  cases 
have  been  performed  to  date,  and  good-to-excellent  agreement  has  been  shown 
between  benchmark  data  and  the  numerical  predictions. 25,26 


4.2  Computer 


For  ail  of  the  computations  completed  in  the  current  study,  a  Silicon  Graphics 
Indigo  was  used.  The  Indigo  had  a  100  mHz  R4000  processor,  and  96  MBytes  of 
RAM,  and  3  GBytes  of  hard  drive  disk  space.  Most  of  the  steady-state  analyses  were 
completed  overnight  on  the  Indigo. 

4.3  Geometry  and  Grid 

The  IFF  configuration  experimentally  tested  earlier  by  UTRC,  and  numerically  by 
CFDRC  in  the  Phase  I  effort  of  this  study,  was  retained  as  the  baseline  configuration 
(Figure  31).  The  model  has  a  0.25  inch  radius  hemispherical  nose  with  a  0.75-inch 
rectangxilar  afterbody,  for  a  total  length  of  1.00  inches,  and  a  height  of  0.5  inches. 
The  total  blockage  was  17%,  with  the  IFF  installed  in  a  3  inch  duct. 

The  computational  domain  started  at  a  station  of  x/H  =  -5,  and  continued  to  x/H  = 
40.  The  time  accurate  grid  consisted  of  four  domains  (Figure  32)  with  the  following 
dimensions: 


zone 

dimensions 

total 

1 

13x53 

689 

2 

30x21 

630 

3 

30x21 

630 

4 

61x77 

4697 

6646 
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Figure  31.  Baseline  IFF  Model  used  for  Current  Study 
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The  last  zone  contained  the  wake  of  the  IFF  and  was  the  most  densely  packed,  with 
37  grid  points  across  the  back  face  of  the  IFF.  For  the  steady-state  solutions,  the  grid 
was  cut  in  half  at  the  centerline,  and  a  symmetry  plane  was  placed  along  the 
centerline. 

All  solutions  used  an  inflow  boundary  condition  set  equal  to  the  experimentally 
measured  flow  values  at  x/H  =  -5.  The  exit  boundary  condition  was  fixed  at  a 
constant  pressure.  The  walls  of  the  IFF  and  the  walls  of  the  tunnel  were  specified  as 
no-slip  boundaries.  Freestream  conditions  were  as  follows: 

U  =24.5  m/s 
P  =  99257  Pa 
T  =389  K 
u/U  =  0.04 
=0.01 

4.4  Numerical  Technique 


Steady  State  Solutions 

Cases  I  and  n  imposed  a  symmetry  plane  at  the  centerline  of  the  geometry, 
preventing  the  asymmetric  vortex  sheddingin  the  wake,  andr^roducing^teady^state 
flow.  The  initial  conditions  were  set  to  the  freestream  conditions,  and  the  solution 
was  converged  4  orders  of  magnitude  for  all  variables.  Convergence  was  typically 
completed  within  1000  iterations,  requiring  6  hours  of  CPU  time. 

After  the  solution  was  completed,  the  variables  of  interest  were  linearly  interpolated 
to  the  experimental  stations  for  comparison  to  the  experimental  data.  The  variables 
used  in  the  data  comparison  were  the  mean  velocity  components  (U,  V),  the  mean 
Reynolds  normal  stresses  (uu,  vv),  and  the  mean  Reynolds  shear  stress  (uv).  For 
the  k-e  turbulence  model,  the  modeled  Reynolds  stress  components  were  extracted 
from  the  solution  using  the  Boussinesq  approximation: 


where: 


^L^  =  0.09pV 


For  Case  n,  the  full  Reynolds  Stress  Model  (RSM)  was  added  to  CFD-ACE.  The  work 
was  completed  by  Dr.Yong  Lai  of  CFDRC,  and  included  several  validation  cases.  A 
summary  of  the  Reynolds  stress  turbulence  model  employed,  and  a  review  of  the 
validation  cases  is  presented  in  Appendix  A. 
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Time-Accurate  Solutions 


Cases  III- VII  involved  solving  the  time-accurate  equations  over  the  full  grid.  Since 
no  symmetry  plane  was  imposed  on  the  centerline,  the  flow  was  free  to  oscillate 
about  the  centerline. 

Solving  the  unsteady  Navier-Stokes  equations  generates  two  contributions  to  the 
turbulence  components.  The  first  is  the  resolved  contribution,  and  is  generated  by 
the  large  scale  fluctuations  in  the  velocity  field  the  CFD  solution  computes.  The 
fluctuations  at  a  given  point  must  be  collected  and  analyzed  with  statistical  methods 
in  a  manner  analogous  to  collected  LDV  data.  The  second  contribution  to  the 
turbulence  comes  from  the  turbulence  model.  The  turbulence  model  contribution 
must  also  be  computed  at  each  time  step  (e.g.  using  the  Boussineq  approximation  for 

the  k-e  model)  and  averaged  over  time.  The  two  contributions  are  added  to  obtain 
the  full  turbulence  contribution. 

Using  turbulence  modeling  with  the  Reynolds  averaged  Navier-Stokes  equations  is 
valid  only  if  the  resolved  mean  flow  fluctuations  and  the  modeled  turbulence 
perturbations  are  not  correlated,  i.e.,  there  exists  a  spectral  gap  in  the  time  scales  of 
the  two  components. 27  If  the  two  components  are  correlated,  the  assumptions  made 
during  the  Reynolds  averaging  process  are  violated,  and  the  more  rigorous  Large 
Eddy  Simulation  (LES)  procedure  must  be  used.  The  LES  procedure  formally 
separates  the  scales  using  a  filtering  technique,  that  restricts  the  turbulence  model  to 
scales  smaller  than  the  resolved  scales. 

The  time  scales  for  the  Reynolds  averaged  approach  can  be  compared  to  determine 
the  validity  of  the  application.  For  the  current  study,  the  resolved  time  scale  can  be 

assumed  to  be  the  inverse  of  the  vortex  shedding  frequency.  For  the  k-e  models,  a 
representative  turbulence  time  scale  (x)  can  be  calculated  by  x  =  k/e. 

For  the  unsteady  cases,  the  solution  procedure  was  started  using  the  same  method  as 
the  steady-state  cases;  the  freestream  conditions  were  used  as  the  initial  conditions, 
and  the  solution  was  run  for  200  iterations  in  the  steady-state  mode.  Since  the 
inflow  conditions  were  not  exactly  symmetric,  the  flow  quickly  began  shedding 
alternating  vortices  in  the  wake.  After  200  iterations,  the  solution  was  switched  into 
the  time-accurate  mode  to  properly  capture  the  vortex  shedding  process. 

In  the  time  accurate  mode,  the  solution  was  progressed  at  time  steps  of  5.0x10-6  sec. 
During  each  time  step,  the  variables  were  converged  at  least  3  orders  of  magnitude. 
The  solution  was  advanced  4500  (~  9  vortex  shedding  cycles)  iterations  in  the  time- 
accurate  mode  to  allow  any  residuals  from  the  initial  start-up  to  exit  the 
downstream  boundarv.  The  solution  was  then  restarted  and  advanced 

J 

approximately  2500  time  steps  to  collect  the  turbulence  statistics.  Since  only  a  few 
vortex  shedding  cycles  were  included  in  the  collection  sample  (~5),  care  was  taken  to 
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ensure  that  the  statistics  were  collected  over  a  whole  period  of  the  vortex  shedding 
cycle.  During  the  statistics  collection  process,  a  sample  was  taken  every  10  time  steps 
and  averaged  with  the  existing  sample  set.  At  the  completion  of  the  run,  the 
appropriate  turbulence  quantities  were  calculated  in  the  same  fashion  as  the 
experimental  values  (see  Appendix  B). 

Case  in  provided  a  solution  to  the  laminar  Navier-Stokes  equations,  using  an 
elevated  laminar  viscosity  (1x10*3  m2/sec)  as  a  simple  model  for  the  effect  of 
tvirbulence.  Such  a  solution  resolves  all  velocity  and  Reynolds  stresses  to  fidelity  of 
the  grid,  without  any  actual  modeling  of  the  turbulence  components. 

Case  IV  employed  the  standard  k-e  turbulence  model  with  the  time-accurate  solver. 
An  evaluation  of  the  time  scales  will  be  presented  with  the  data,  to  assess  the 
_ presence-of  a  spectral  gap-between  the  resolved  and  modeled  tmbulence  scales. _ 

Case  V  used  a  relatively  new  formulation  of  the  k-e  turbulence  model  derived  from 
the  Renormahzed  Group  theory  (RNG).  The  RNG  k-e  model  is  functionally  very 
similar  to  the  standard  k-e  model.  The  major  difference  is  that  the  ad-hoc  constants 
used  in  the  standard  model  are  now  derived  from  the  RNG  theory.  The  differences 
in  the  constants  are  as  follows: 

Standard  k-e  model 

Cf,  =  0.09,  Cel  =  2-44,  Cez  =  I  -92 
(Tk  =  l-0,  <T£=1.3 

RNG  k-e  model 


C  =0.085, 
Cel  =  1-68, 


'el 


1.42- 


’(1  - 

l+p7]^ 


Gk  =  0.7179,  Oe  =  0.7179 


=438  |3  =  0.015 

O 


The  RNG  k-e  model  is  derived  by  making  an  expansion  about  an  equilibrium  state 
with  known  Gaussian  statistics,  with  the  effects  of  mean  strains  being  represented  by 
a  random  force.  The  bands  of  the  smallest  scales  are  systematically  removed  and  the 
remaining  space  is  rescaled.  28  This  process  is  continued  until  only  the  resolvable 
scales  remain.  Speziale  and  Thangam29  have  shown  excellent  resvdts  for  a  backward 
facing  step  using  the  RNG  k-e  model.  While  the  RNG  model  can  easily  be 
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integrated  down  to  the  wall  without  the  use  of  wall  functions,  the  method  used  in 
the  current  study  still  employed  the  standard  wall  functions. 

Cases  VI  and  Vn  were  also  unsteady  and  used  the  RNG  k-e  model,  but  solved  the 
3-D  equations  to  provide  an  extra  degree-of-freedom  to  the  resolved  turbulence 
component.  Case  VI  used  a  relatively  small  spanwise  domain  (~  IH)  while  Case  VII 
covered  a  larger  spanwise  domain  (~  6H). 

4.5  Numerical  Results 


Steady  State 

The  results  for  the  two  steady  state  cases  are  presented  in  Figures  33-37,  and 
compared  to  the  experimental  data.  The  results  are  shown  for  the  near  wake  region 
only,  from  x/H  =0.5  to  x/H  =1.2.  The  predicted  mean  velocities  for  Case  I  indicate 
that  the  numerical  recirculation  bubble  for  both  cases  extends  to  x/H=2,  much 
further  downstream  than  the  experimentally  determined  bubble  position.  The 
numerical  bubble  length  of  two  base  heights,  or  four  half  bluff  body  heights,  is 
similar  to  that  of  a  backward-facing  step.  This  is  to  be  expected  since  the  model 
resembles  a  backward-facing  step  with  a  slip  wall  at  the  centerline.  The  predicted 
Reynolds  stresses  are  also  comparable  to  those  measured  in  backward-facing  step 
flows.  However,  the  experimentally  measured  Reynolds  stresses  for  the  IFF  were 
found  to  be  ~5  times  higher  than  those  found  in  backward-facing  step  flows. 
Additionally,  the  turbulence  values  were  clearly  nonisotropic  for  the  IFF  flow,  and 
the  k-e  model  predicted  isotropic  turbulence. 

The  Reynolds  stress  model,  while  not  restricted  to  isotropic  turbulence,  does  not 
show  a  marked  improvement  over  the  standard  k-e  model.  The  RSM  solves  four 
more  nonlinear  partial  differential  equations  than  the  k-e  model  and  requires 
considerably  more  time  to  converge.  In  the  case  of  the  flameholder,  the  RSM  model 
is  not  worth  the  increased  computational  effort.  Since  the  model  does  not  improve 
the  solution  by  any  significant  amount,  the  error  is  probably  not  due  to  the  nature  of 
isotropic  turbulence. 
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Figure  33- 


Normalized  Mean  Axial  Velocity  and  V 
(O  data  - RSM  ■  •  *  k-e) 
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Figure  35.  Normalized  uu  Turbulent  Stress  and  Stress  Gradient  Distributions 
(O  data  ■  -  *  k-e) 
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Figure  36.  Normalized  vv  Turbulent  Stress  and  Stress  Gradient  Distributions 
(o  data  ■  *  -  k-e) 


58 


3  fTTTrprrTprrrrp- 

L 

h  : 

2H 
r 

1  H 


4- 


T 

I 


t  0^ 

r 

r 

•1 

F 

C 

-2  b 


t 

4- 

T 

T 

T 


JL  : 


T 


.3  * 


+ 

I 


t  -r 


j. 

+ 


± 


t  J 


T 


Juiiiliulr 


-1.0  0.0  1.0 


(H/l/2)auv/9x 


Figure  37.  Normalized  uv  Turbulent  Stress  and  Stress  Gradient  Distributions 
(O  data  ■  -  -  k-e) 
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Time  Accurate 


The  time  accurate  solutions  were  initiated  in  an  attempt  to  improve  upon  the 
steady-state  computational  results.  Since  the  experimental  results  indicated  that  the 
IFF  was  shedding  alternating  large-scale  vortices  from  the  top  and  bottom  corners  of 
the  IFF,  significant  turbulence  must  be  generated  from  the  vortices.  Such 
turbulence  would  be  beyond  that  generated  by  the  backward-facing  step.  Therefore, 
in  order  to  properly  predict  the  higher  turbulence,  the  unsteady  nature  of  the 
flowfield  must  be  taken  into  account.  It  is  not  reasonable  to  believe  that  a  single 
turbulence  model  could  properly  predict  a  backward  facing  step  and  a  bluff  body, 
since  the  only  geometrical  difference  is  the  boundary  layer  at  the  symmetry  plane.  A 
time-accurate  simulation  is  necessary  to  properly  capture  the  increased  hirbulence. 

The  three  time-accurate  calculations  (Cases  HI,  IV,  and  V)  were  able  to  predict  the 
vortex  shedding  behind  the  IFF,  and  a  sample  result  showing  the  instantaneous 
streamlines  is  presented  in  Figure  38.  The  structure  of  the  vortex  shedding  is 
evident  as  the  vortex  from  the  lower  surface  has  been  shed  downstream  and  the 
upper  surface  vortex  is  forming.  The  vortex  pattern  continues  downstream,  but  as 
the  vortex  convective  velocity  accelerates  to  match  the  freestream  velocity,  the 
presence  of  the  vortices  in  the  streamlines  is  disguised. 


Figure  38.  Streamlines  Show  Unsteady  Vortex  Motion  Behind  Example  Bluff  Body 
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Case  ni  used  a  constant  viscosity  of  1x10-3  m2/sec,  and  predicted  a  Strouhal  number 
of  0.271,  15%  higher  than  the  experimental  value  of  0.243.  The  mean  and 
fluctuating  velocity  components  are  compared  in  Figures  5.9-5.13.  The  velocity  data 
from  the  laminar  unsteady  solution  showed  much  better  agreement  with  the  data 
than  the  steady  state  results.  The  predicted  recirculation  bubble  length  of  1.08  H  was 
much  closer  to  the  experimental  value  of  0.9  H.  In  addition,  the  qualitative 
agreement  of  the  turbulence  profiles  indicated  that  the  vortex  shedding  is  the 
dominant  contributor  to  the  turbulence  in  the  experiment.  A  close  examination  of 
the  data  revealed  that  the  strength  of  the  shed  vortices  was  probably  too  strong.  To 
check  the  effect  of  the  laminar  viscosity  on  the  solution,  the  run  was  repeated  for 
two  different  values  of  viscosity,  1x10-2  m2/sec  and  1x10-5  in2/sec.  At  the  higher 
value  of  viscosity,  the  vortex  shedding  was  damped  and  the  solution  returned  to  the 
steady-state  results.  At  the  lower  viscosity,  the  vortex  shedding  frequency  was 
increased  slightly,  but  the  strength  of  the  vortices  increased,  significantly  degrading 
the  data  comparison.  Therefore,  it  was  concluded  that  more  accurate  turbulence 
modeling  was  required  to  properly  predict  the  flowfield. 

Case  IV  addressed  the  issue  of  turbulence  modeling  by  employing  the  standard  k-e 
turbulence  model  with  the  standard  wall  functions.  The  first  attempt  at  using  the  k- 
e  turbulence  model  with  the  unsteady  flow  yielded  a  steady-state  solution.  Other 
researchers30  have  found  similar  results  and  have  blamed  the  numerical 
formulation  of  the  turbulence  model.  However,  later  studies  found  the  damped 
behavior  to  be  a  function  of  the  differencing  scheme  employed  instead  of  the 
turbulence  model.3i  The  first  attempt  in  the  current  study  employed  first  order 
differencing  in  space.  When  central  differencing  was  used,  the  vortex  shedding  was 
maintained.  The  k-e  model  predicted  the  closest  value  for  the  Strouhal  number  at 
0.240  (~1%  low).  The.  mean  velocities  and  Reynolds  stresses  are  shown  in  Fibres 
39-43.  Examination  of  the  mean  axial  velocities  indicates  that  the  k-e  solution  is  not 
much  different  from  the  laminar  case,  with  the  length  of  the  recirculation  zone 
about  the  same.  However,  the  transverse  velocity  does  show  some  differences. 
Very  near  the  IFF,  the  k-e  model  shows  a  better  match  to  the  data,  with  the 
inflection  point  at  the  centerline  somewhat  reduced.  But  further  downstream,  the 
laminar  case  maintains  the  higher  transverse  velocities  reflected  by  the 
experimental  data.  This  is  due  to  the  higher  turbulent  viscosity  predicted  in  the 
wake  by  the  k-e  model  (0.004  vs  0.001  for  the  laminar  case).  The  normal  axial  stress 
was  also  similar  to  the  laminar  case  very  near  the  IFF.  But  the  axial  trends  were  not 
correct  for  either  the  k-e  or  laminar  cases,  showing  increasing  peak  stresses  with 
axial  position,  while  the  experiment  showed  a  decreasing  peak  stress. 
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Figure  39.  Normalized  Mean  Axial  Velocity  and  Velocity  Gradient  Distributions 
(O  data  - RNG-2D  '■■laminar - k-e) 
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Figure  41.  Normalized  uu  Turbulent  Stress  and  Stress  Gradient  Distributions 
(O  data  - RNG-2D  ■  •  •  laminar - k-e) 


Figure  43.  Normalized  uv  Turbulent  Stress  and  Stress  Gradient  Distributions 
(O  data  - RNG-2D  -••■laminar - k-e) 
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In  an  effort  to  improve  upon  the  predictions  of  the  k-e  turbulence  model,  a  different 
derivation  of  the  k-e  model  based  on  Renormalized  Group  theory  (RNG),  was  tested 
(Case  V).  The  RNG  k-e  model  is  functionally  very  similar  to  the  standard  k-e 
model,  except  that  the  constants  are  derived  from  the  theory  instead  of  being 
empirically  correlated  with  experimental  data.  In  addition,  one  of  the  constants, 
Cgi,  is  no  longer  constant,  but  a  function  of  local  strain  rates. 

The  results  for  Case  V  are  also  presented  in  Figures  39-43.  The  RNG  k-e  results  were 
the  best  of  the  group,  yielding  very  good  results  for  the  axial  mean  velocities. 
Although  the  transverse  velocities  compared  well  at  the  initial  axial  stations,  the 
comparison  degraded  with  increasing  axial  position,  indicating  that  the  vortex 
structures  were  dissipating  too  quickly  in  the  numerical  simulation.  Both  of  the 
normal  stresses  were  overpredicted  by  the  RNG  k-e  model,  but  the  axial  trends  were 
correct.  The  same  can  be  said  for  the  shear  stress;  the  numerical  values  were 
slightly  overpredicted,  but  the  trends  were  correct.  In  both  cases,  the  fact  that  the 
trends  were  correct  is  significant,  since  with  the  laminar  and  standard  k-e  cases,  the 
trends  were,  in  general,  not  correct.  The  RNG  k-e  model  overpredicted  the  Strouhal 
number  by  4%  (0.253). 

After  the  completion  of  the  time-accurate  solutions,  an  analysis  of  the  relative  time 
scales  was  possible.  In  all  cases,  the  vortex  shedding  frequency  was  near  450  Hz, 
yielding  a  mean-flow  time  scale  of  x  =  2x10-3  sec.  The  standard  k-e  model  produced  a 
maximum  time  scale  in  the  wake  of  x  =  2x10-3  sec.  The  fact  that  the  time  scales  are 
the  same  indicates  that  a  formal  Large  Eddy  Simulation  (LES)  is  required.  However, 
the  region  within  the  recirculation  zone  had  a  much  lower  k-e  time  scale  (2x10-4 
sec),  indicating  that  the  near  field  is  probably  accurate.  The  RNG  k-e  model 
generated  similar  time  scales. 

The  results  presented  in  the  current  study  are  not  unlike  that  of  other  researchers. 
Martensson,  et  al.32^  obtained  similar  results  using  a  Large  Eddy  Simulation  (LES) 
technique,  and  found  that  the  transverse  normal  stress  was  reduced  significantly 
when  simulated  with  a  3-D  analysis.  Two  3-D  analyses  were  completed  as  part  of  the 
present  study,  both  using  the  RNG  k-e  turbulence  model.  The  first  analysis 
involved  a  narrow  slice  of  the  tunnel,  spanning  only  one  base  height  across  the 
tunnel.  The  second  solution  utilized  a  larger  spanwise  section,  covering  6  base 
heights  across  the  timnel. 

For  the  purpose  of  the  3-D  study,  the  grid  was  reduced  in  density  to  allow  better 
resolution  in  the  spanwise  direction  with  the  available  computer  resources.  The 
new  grid  had  the  following  dimensions  and  is  shown  in  Figure  44. 
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Figiire  44.  Coarse  Grid  used  for  3-D  Unsteady  Analysis 


zone 

dimensions 

total 

1 

9x41x25 

9225 

2 

20  X  15  X  25 

7500 

3 

20  X 15  X  25 

7500 

4 

41  X  65  X  25 

66625 

90850 

To  test  the  effects  of  the  new  grid,  a  2-D  case  was  repeated  for  the  coarse  mesh  and 
the  results  compared  to  the  those  from  the  original  grid.  At  the  near  field  stations, 
the  results  were  virtually  identical.  But  by  2  base  heights  downstream  of  the  IFF,  the 
comparison  began  to  degrade  at  the  centerline,  producing  a  5%  error  in  U.  Based  on 
these  results,  the  new  coarse  grid  was  considered  sufficient  to  study  the  effects  of  the 
3-D  analysis  on  the  near  field  solution  only. 

The  first  3-D  solution  (Case  VI)  used  periodic  boundary  conditions  on  the  spanwise 
boundaries,  and  covered  1  base  height  with  25  planes  in  the  spanwise  direction.  The 


resulting  time-averaged  flowfield  closely  resembled  that  of  the  symmetric  steady- 
state  solution,  with  a  long  recirculation  zone.  Arnal  and  Freidrich33  have  reported 
similar  difficulties  with  a  backward  facing  step,  indicating  that  spanwise  extensions 
less  than  4  step  heights  can  artificially  elongate  the  separation  bubble.  Why  this 
happens  is  not  clear. 

As  a  result,  a  new  solution  (Case  VII)  was  completed  covering  a  larger  spanwise 
domain.  The  new  grid  covered  6  base  heights,  ^so  using  25  spanwise  planes  and 
periodic  spanwise  boundary  conditions.  One  major  difference  was  noted  in  the 
behavior  of  the  2-D  and  3-D  solutions.  The  2-D  solutions  established  a  steady 
periodic  motion  very  quickly,  typically  within  0.005  seconds.  The  3-D  solutions 
went  through  several  different  modes  before  settling  down  to  what  appeared  to  be  a 
higher  frequency  mode  (-500  Hz)  superimposed  on  a  lower  frequency  mode  (-100 
Hz).  Therefore  the  3-D  solutions  had  to  be  run  over  a  much  longer  period  of  time  to 
establish  the  periodic  motion.  In  addition,  the  amount  of  data  collected  to  calculate 
accurate  statistics  increased  dramatically,  since  the  lowest  frequency  dictated  the 
sampling  window.  The  results  from  Case  Vn  are  shown  in  Figures  45-49.  There  is  a 
slight  improvement  in  the  mean  velocities  at  all  stations  as  compared  to  the  2-D 
results.  The  predicted  recirculation  zone  length  matches  the  experimental  results  to 
within  the  fidelity  of  the  experiments  (x/H  =  0.9).  However,  the  predicted  Strouhal 
number  is  now  7%  below  the  experimental  value  0.226  vs.  0.243).  The  Reynolds 
stresses  were  reduced  for  all  values,  and  compared  very  well  with  the  experimental 
values.  The  exception  was  the  transverse  normal  stress,  which  was  reduced  from 
the  2-D  results,  but  was  reduced  well  below  the  experimental  results. 

The  effect  of  the  added  degree  of  freedom  from  the  third  axis  of  motion  is  several 
fold.  First,  the  three  dimensional  vortex  filament  shed  from  the  IFF  is  now  free  to 
distort  in  the  spanwise  direction.  The  low  frequency  made  in  the  3-D  solution  is 
evidence  of  such  distortions.  As  a  result  of  the  spanwise  distortions,  the  shed 
vortices  will  make  a  contribution  to  the  spanwise  stresses,  dependant  upon  the 
magnitude  of  the  spanwise  distortion.  Any  contribution  to  the  spanwise  stresses 
will  result  in  a  likewise  decrease  to  the  axial  and  transverse  stresses,  dependant 
upon  the  local  direction  of  the  distortion.  Such  an  idea  was  put  forth  by  Perry  and 
Steiner34,  and  is  supported  by  the  present  results. 

The  second  effect  of  the  vortex  filament  distortion  is  to  create  a  vortex  stretching 
effect.  As  the  vortex  filament  is  distorted,  the  total  length  is  increased.  As  a  result, 
the  vortex  velocity  increases  to  maintain  the  same  angular  momentum.  Such 
vortex  stretching  could  explain  the  consistent  lack  of  agreement  between  the  mean 
transverse  velocity  components  moving  downstream  of  the  IFF. 
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Figure  45.  Normalized  Mean  Axial  Velocity  and  Velocity  Gradient  Distributions 
(O  data  - RNG-2D  - RNG-3D) 


Figure  48.  Normalized  uu  Turbulent  Stress  and  Stress  Gradient  Distributions 
(O  data  - RNG-2D  - RNG-3D) 
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5.0  ADVANCED  MODEL  NUMERICAL  TESTS 


Running  concurrently  with  the  baseline  tests  was  an  effort  to  find  a  better  design  for 
the  IFF.  The  goal  was  to  create  a  flameholder  with  increased  fuel /air  mixing, 
without  paying  a  high  penalty  in  pressure  loss.  On  some  of  the  advanced  models, 
simple  combustion  models  were  analyzed  to  establish  the  improvement  to  the 
combustion  efficiency  based  on  increased  mixing.  On  all  of  the  models,  the  cold 
pressure  loss  was  calculated  for  comparison  to  the  baseline. 

A  total  of  27  different  models  were  numerically  tested  in  this  phase  of  the  study.  A 
variety  of  different  approaches  were  taken  in  an  attempt  to  generate  the  large  scale 
axial  vortices.  On  the  different  models,  58  production  CFD  runs  were  completed, 
and  countless  smaller  test  cases  were  also  completed.  In  the  following  sections,  the 
computational  methods  used  for  the  study  will  be  reviewed.  A  summary  of  all  the 
models  tested  will  then  be  presented  followed  by  a  comparison  of  the  data  obtained 
from  the  numerical  tests. 

5.1  Computational  Method 

The  same  CFD  code  was  used  for  the  advanced  model  testing,  that  was  used  for  the 
baseline  isothermal  tests  (CFD-ACE,  see  Section  3.1).  The  major  difference  is  that 
most  of  the  solutions  for  the  advanced  models  were  completed  using  the  steady 
state  solver.  While  ignoring  the  unsteady  contributions  may  have  caused  a  small 
loss  of  accuracy,  the  trends  were  still  evident  between  the  different  models. 

Instantaneous  reactions  were  used  for  all  combustion  calculations,  based  on  the 
following  stoichiometric  equation: 

C3H8  +  5  O2  ->  3  CO2  +  4  H2O 

It  was  felt  that  the  simple  model  would  identify  the  proper  trends  between  the 
different  models.  For  the  rigorous  combustion  tests  (presented  in  Section  8),  more 
advanced  combustion  models  were  used. 

The  combustion  efficiency,  as  used  in  the  current  study,  is  simply  a  fimction  of  the 
mixing  occurring  between  the  fuel  and  the  air.  Combustion  efficiency  is  defined  as: 

heat  released  ^  fuel  tmrnea 
Icomb  available  in  fuel  total  fuel 

For  a  mixing-controlled  combustion  model,  combustion  efficiency  will  always  be 
100%  in  any  grid  cell  that  is  fuel  lean.  Efficiencies  less  than  100%  will  only  occur  in 
grid  cells  that  are  fuel  rich,  and  represent  the  presence  of  unbumed  fuel.  In  the 
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current  program,  combustion  efficiency  was  postprocessed  by  mass  weighting  and 
summing  the  combustion  efficiencies  in  each  cell  in  a  specified  axial  plane. 

The  pressure  losses  for  the  flameholders  models  must  be  computed  based  on  cold 
flow  results.  The  heat  release  of  combustion  adds  to  the  pressure  loss,  and 
complicates  the  comparison  of  different  models.  If  one  model  has  better 
combustion  efficiency  &an  another,  it  will  also  have  more  pressure  loss  due  to  heat 
release.  Only  cold  flow  pressure  losses  can  be  directly  compared.  Since  the  presence 
of  the  crossflow  jet  may  also  contribute  to  the  pressure  loss,  all  pressure  loss 
calculations  were  completed  with  solutions  that  included  the  fuel  injection,  but  did 
not  allow  combustion. 

The  total  pressure  loss  associated  with  a  bluff  body  can  be  computed  in  one  of  two 
ways.  First,  the  mass  weighted  total  pressures  in  each  cell  can  be  integrated  across  an 
axial  plane  to  give  a  direct  calculation  of  the  total  pressxire.  Or,  the  static  pressure 
can  simply  be  evaluated  at  the  side  wall  of  the  duct  downstream  of  the  bluff  body, 
and  represents  a  good  approximation  to  the  total  pressure  loss  when  compared  with 
the  static  pressure  upstream.  When  using  the  second  method,  the  static  pressure 
must  be  evaluated  far  enough  downstream,  since  the  static  pressure  will  partially 
recover  in  the  wake  of  the  bluff  body. 

The  numerical  solutions  were  completed  using  first  order  differencing  in  space  on 
all  variables.  To  model  turbulence  effects,  the  standard  k-e  turbulence  model  was 
used  along  with  the  standard  wall  fimctions.  All  of  the  solutions  were  considered 
converged  when  all  residuals  were  reduced  by  at  least  four  orders  of  magnitude. 

5.2  Inlet  Conditions 

The  freestream  inlet  conditions  were: 

U  =  66  m/s 
M  =  0.2 

P  =  20700  N/m2 
T  =  272  K 
u/U  =  0.10 

Pt  =  0.01 

working  gas:  air 

Uniform  inlet  conditions  were  used  in  all  cases,  and  the  inlet  plane  was 
approximately  20  base  heights  upstream  of  the  IFF  to  allow  a  boundary  layer  to 
develop. 


In  the  cases  that  included  fuel  injection,  the  following  jet  conditions  were  enforced: 
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V  =  114  m/s 
M  =  0.42 
P  =  20700  N/m2 
T  =  272 
u/U  =  0.05 
=  0.01 

working  gas:  propane 

Unless  otherwise  noted,  the  jet  conditions  were  set  up  to  produce  an  equivalence 
ratio  of  0.5  and  a  jet-to-freestream  momentum  flux  ratio  (J)  of  4.4. 

53  Advanced  Flameholder  Models 

53.1  Baseline  Model 

Before  any  advanced  models  were  tested,  a  new  baseline  had  to  be  established.  The 
new  baseline  was  different  from  the  baseline  tests  presented  in  Section  4,  due  to 
several  differences  in  the  cases.  First,  the  new  baseline  was  three-dimensional  since 
it  included  the  fuel  jets  in  the  solution.  Second,  the  inlet  plane  was  set  to  a  uiuform 
velocity.  Finally,  the  new  baseline  also  needed  to  include  combustion  to  determine 
the  baseline  combustion  efficiency  curve. 

The  grid  for  the  new  baseline  case  contained  36,465  grid  points  (55x39x17)  and  is 
shown  in  Figure  50.  The  grid  covers  only  the  top  half  of  the  domain,  and  spans 
from  the  midpoint  of  an  injector  to  the  midplane  between  injectors.  Symmetry 
planes  were  used  on  both  spanwise  planes  and  at  the  tuimel  centerline  plane. 
Viscous  boundaries  were  established  on  the  top  hmnel  wall,  as  well  as  on  the  walls 
of  the  flameholder.  The  physical  dimensions  of  the  flameholder  were  the  same  as 
those  in  the  isothermal  tests  (Section  4)  .  The  fuel  injector  measured  0.115  inches  in 
diameter,  and  was  modeled  with  a  4x4  round  grid  pattern.  Injectors  were  vertically 
aligned  and  were  spaced  0.5  inches  apart.  The  baseline  geometric  blockage  was 
16.7%.  The  grid  is  recognized  as  being  coarse,  and  a  finer  grid  was  run  to  ensure  the 
coarse  grid  was  sufficient  to  capture  qualitative  features  of  the  flowfield.  The  fine 
grid  solution  compared  well  with  the  coarse  grid  solution;  in  order  to  reduce  run 
times,  the  coarse  grid  was  selected  for  parametric  studies. 
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5.3.2  Advanced  Model  ADV2 

The  first  advanced  model  tested  was  called  ADV2.  The  name  actually  applied  to  an 
entire  model  group  with  several  parametrics  analyzed  within  the  group.  The  first 
model  tested  was  ADV2  aO,  and  is  shown  in  Figure  51.  The  total  flameholder 
length  was  1.34  inches  with  a  maximum  height  of  ±0.75  inches.  This  model 
followed  closely  the  Split-Notched  IFF  model  from  the  Phase  I  study.  It  used  square 
fingers  that  protruded  into  the  flow  from  the  trailing  edge.  The  base  area  of  each 
finger  was  equal  to  the  original  base  area  of  the  baseline  model,  therefore  the  total 
base  area  at  the  trailing  edge  (20.8%)  was  only  increased  by  the  extra  material 
forming  a  web  between  the  fingers.  However,  the  geometric  blockage  created  by  the 
possible  flow  separating  off  of  the  ramps,  was  much  larger  at  35.4%.  The  CFD 
analysis  was  not  set  up  to  predict  separation,  due  to  the  grid  density  and  turbulence 
modeling,  although  the  code  did  show  some  separation  within  the  channel  regions. 
The  ramps  were  placed  at  45  degrees  on  both  top  and  bottom,  and  were  expected  to 
encounter  some  flow  separation. 

The  second  model  from  the  ADV2  group  was  called  a  10,  and  was  generated  by 
slanting  the  web  between  the  fingers  at  a  10  degree  angle  (Figure  52).  The  intent  was 
to  increase  the  axial  vorticity  by  providing  a  sharper  edge  to  generate  the  vortex. 
The  most  important  effect  of  angling  the  web,  was  the  increase  in  base  area  due  to 
the  stretching  of  the  finger  base  areas.  The  new  base  area  was  23.8%,  with  a 
geometric  blockage  of  38.4%.  All  other  aspects  of  the  flameholder  were  the  same  as 
the  aO  model. 

The  same  idea  of  angling  the  web  between  the  fingers  was  carried  further  with  the 
next  two  models,  a20  and  a30  (Figures  53  and  54).  The  a20  model  had  a  base 
area /geometric  blockage  of  26.9%/41.5%  while  the  (x30  model  had  30.5%/ 45.0%. 

The  results  from  the  combustion  tests  are  shown  in  Figure  55.  From  the 
combustion  efficiency  curves,  each  advanced  model  can  be  seen  to  produce  faster 
mixing  than  the  baseline.  This  is  due  to  the  formation  of  a  large  axial  vortex 
downstream  of  the  flameholder,  as  can  be  seen  in  Figure  56.  Of  the  four  models 
tested,  the  a20  model  showed  slightly  better  overall  performance.  Therefore,  it  was 
chosen  for  further  parametric  studies. 
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Figure  51.  3  D  Model  of  ADV2  aO,  and  Definition  of  a, 


Figure  52.  3-D  Model  of  ADV2  alO 


Figure  54.  3-D  Model  of  ADV2  (x30 
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Figure  55.  Combustion  Efficiency  for  ADV2  a  Models  Shows  a20  Gives  Superior  Performance 
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From  Figure  57,  the  presence  of  a  strong  fuel  core  can  be  seen  above  the  ramp.  The 
fuel  is  from  the  injector  directly  in  front  of  the  ramp.  Possible  performance  gains 
were  recognized  if  the  fuel  core  could  be  better  distributed.  The  first  trial  at 
redistributing  the  fuel  core  was  to  change  the  angle  on  the  top  of  the  ramp.  Two 
different  angles  were  tested:  P20  and  [3-20  shown  in  Figures  58  and  59.  The  intent 
was  to  create  a  secondary  vortex  located  near  the  top  of  the  ramp.  However,  no 
vortex  was  formed  and  the  fuel  distribution  was  not  significantly  changed.  The  new 
combustion  efficiency  curves  are  shown  in  Figure  60,  and  show  that  the  P20  model 
does  not  change  the  performance  much,  even  though  the  blockage  is  considerably 
higher  (32.8%/ 47.4%).  The  P-20  model  shows  reduced  performance  with  blockages 
of  (19.8%/34.4%). 

The  next  attempt  to  redistribute  the  fuel  core  was  to  eliminate  the  fuel  jet  directly  in 
front  of  each  ramp,  thereby  eliminating  the  fuel  core  (model  JETMOD).  The  fuel 
lost  to  the  removed  injectors  is  replaced  by  additional  fuel  introduced  through  the 
remaining  injectors.  To  keep  the  momentum  flux  ratio  the  same,  the  increased  fuel 
flow  was  introduced  by  increasing  the  size  of  the  remaining  injectors  to  a  diameter 
of  0.133  inches.  The  trial  was  modeled  on  the  ADV2  a20  model  only.  The  fuel 
distribution  is  shown  in  Figure  61  for  the  two  cases.  Although  the  fuel  core  has 
been  removed,  the  remaining  fuel  is  slightly  more  concentrated  in  the  center  of  the 
flow.  The  resulting  combustion  efficiency  curves  are  shown  in  Figure  62,  where  the 
modification  is  seen  to  have  little  if  not  a  detrimental  effect.  The  ramps  clearly 
serve  two  purposes.  Not  only  do  they  provide  for  the  strong  axial  vortidty,  but  they 
also  direct  the  fuel  to  the  outer  regions  of  the  flow  above  and  below  the 
flameholder. 

Of  all  the  models  tested  in  the  preceding  series,  ADV2  a20  showed  the  best 
combustion  performance.  However,  the  pressure  loss  associated  with  the  mixing 
enhancement  must  also  be  considered.  Figure  63  shows  the  cold  pressure  loss  of  all 
the  models  tested.  Clearly  the  pressure  losses  were  much  higher  than  the  baseline 
IFF.  Even  a  double  base  IFF  (two  IFF's  placed  vertically  in  tandem)  had  a  much 
lower  pressure  loss,  while  providing  the  most  flame  stability,  based  on  the  base 
area. 
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Figure  57.  Strong  Fuel  Core  Forms  Above  Ramp;  ADV2  a20 
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Figure  60.  Combustion  Efficiency  for  ADV2  p  Models 


Figure  61.  Model  JETMOD  Reduces  Presence  of  Fuel  Core 
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Figure  62.  Model  JETMOD  does  not  Improve  Combustion  Efficiency 


single  IFF 
double  IFF 
adv2a0 

adv2a0  34  ramp 
adv2a0  26  ramp 
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Figure  63a.  Pressure  Loss  for  All  Models  Analyzed 
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adv5 

adv6 
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Figure  63b.  Pressure  Loss  for  All  Models  Analyzed 


In  an  effort  to  reduce  the  pressure  loss,  a  parametric  study  was  completed  on  the 
effect  of  reducing  the  ramp  angle.  To  keep  the  base  area  constant,  the  tail  of  the 
flameholder  was  axially  stretched  until  the  desired  angle  was  attained.  Four 
different  tail  lengths  were  tested,  using  the  aO  model  as  a  baseline  (see  Figure  64).  In 
this  study,  the  maximum  crossflow  velocity  at  a  given  station  during  nonreacting 
tests,  iirstead  of  the  combustion  efficiency,  was  used  to  judge  the  mixing  ability  of 
the  different  designs.  The  results  are  shown  in  Figure  63,  where  the  pressure  loss  is 
seen  to  be  directly  related  to  the  ramp  angle;  the  lower  the  angle,  the  lower  the 
pressure  loss.  However,  the  crossflow  velocity  does  not  follow  the  same  trend.  The 
34  degree  ramp  provides  the  highest  crossflow  velocity,  but  still  a  lower  pressure 
loss  than  the  45  degree  baseline.  As  the  ramp  angle  decreases  further,  the  crossflow 
velocity  declines.  Based  on  these  comparisons,  the  best  choice  based  on  both  criteria 
would  be  the  26  degree  ramp,  having  a  substantial  crossflow  component  and  a 
relatively  low  pressure  loss. 

As  a  final  check  with  the  ADV2  series,  the  aO  model  was  tested  with  and  without  a 
web  between  the  fingers.  The  results  showed  a  small  decrease  in  pressure  loss  (0.4%) 
and  virtually  no  difference  in  the  crossflow  profile.  Therefore,  there  was  a  slight 
performance  gain  by  not  providing  a  web  between  the  fingers.  Removing  the  web 
also  removed  an  area  that  is  difficult  to  cool,  due  to  the  relative  thickness  of  the 
web.  However,  removing  the  web  removes  any  connection  between  the  base  faces 
of  individual  fingers.  Such  connectivity  may  be  necessary  for  proper  ignition 
propagation  across  the  flameholder. 

5.3.3  Advanced  Model  ADV3 

The  next  advanced  model,  ADV3,  took  a  slightly  different  approach.  The  previous 
models  incorporated  ,  large  base  areas  at  the  end  of  each  finger  with  a  thin  web 
connecting  the  fingers.  ADV3  took  the  opposite  approach  by  providing  a  think  web, 
with  thin  connecting  regions  at  the  ends  of  each  finger  (see  Figure  65).  By 
removing  much  of  the  material  from  the  end  of  the  fingers,  the  fingers  did  not 
protrude  as  far  into  the  freestream,  and  thus  reduced  the  frontal  area  of  the 
flameholder.  The  side  angle  of  a20  was  retained  from  the  previous  tests,  and  the 
inward  ramps  were  kept  at  45  degrees,  and  had  the  same  possible  separation 
problems  as  the  ADV2  models. 

The  results  from  ADV3  are  shown  in  Figures  66  and  63.  The  pressure  losses  were 
greatly  reduced  compared  to  the  other  models,  due  to  the  reduced  blockage. 
However,  the  combustion  efficiency  was  considerably  less  than  the  ADV2  a20 
model.  In  addition,  the  characteristic  dimension  of  the  flameholder  face  was  only 
0.25  inches  -  half  that  of  ADV2  -  and  therefore  would  provide  less  flame  stability. 
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Combustion  Efficiency  for  ADV3  Models 


A  second  variant  of  ADV3  was  tested,  and  was  generated  by  extending  the  trailing 
edge  downstream  another  0.375  inches,  while  holding  the  inner  ramp  angle 
constant  at  45  degrees  (Figure  67).  The  outer  web  thickness  was  also  held  constant  at 
0.125  inches,  providing  a  30  degree  outward  facing  ramp.  The  results  are  also  shown 
in  Figures  66  and  63.  The  combustion  efficiency  was  very  close  to  that  of  the  ADV2 
a20  model.  But  the  geometric  blockage  is  comparable  to  that  of  ADV2  a20,  and  the 
pressure  loss  reflects  that  fact,  matching  the  ADV2  a20  results. 

5.3.4  Advanced  Model  ADV4 

The  advanced  model  ADV4  was  the  first  of  two  exotic  models  intended  to  explore 
different  ways  of  using  the  momentum  of  the  fuel  to  help  create  the  axial  vortices. 
ADV4  is  shown  in  Figure  68  and  was  designed  from  the  extended  version  of  ADV3, 
by  joining  the  tips  of  the  adjacent  fingers  together  with  a  fuel  bar.  The  fuel  was  then 
injected  through  the  fuel  bar  instead  of  through  the  injectors  located  in  the  mid¬ 
section  of  the  IFF.  The  goal  was  to  use  the  fuel  to  enhance  attached  flow  passing 
through  the  channels  between  the  fingers  (see  Figure  69).  With  the  attached  flow, 
the  axial  vortex  would  be  strengthened.  The  results  shown  in  Figures  70  and  63 
show  that  the  design  performed  poorly  in  all  aspects.  Having  the  highest  geometric 
blockage  (47%),  the  pressure  losses  are  high.  In  addition,  the  fuel  is  injected  from 
below  the  centerline,  degrading  fuel  penetration.  As  a  result,  the  mixing 
performance,  and  combustion  efficiency,  is  very  low. 


5.3.5  Advanced  Model  ADV5 

The  next  advanced  model,  ADV5,  was  also  an  exotic  approach.  ADV5  relied  on  the 
principles  of  the  Coanda  effect  and  slotted  flap  technology  to  increase  the  mixing 
performance.  The  Coanda  effect  has  been  demonstrated  to  reduce  the  drag  of  bluff 
bodies. 35  Shown  in  Figure  71,  ADV5  injects  the  fuel  tangentially  with  the 
freestream,  using  the  Coanda  effect  to  pull  the  flow  around  the  curved  surface  of  the 
finger.  The  remaining  surface  of  the  finger  was  energized  with  the  high 
momentum  flow  coming  through  the  slotted  flap  located  below  the  fuel  injector. 

A  two-dimensional  model  was  first  tested  to  verify  the  validity  of  the  design.  Figure 
72  shows  the  velocity  vectors  that  resulted  from  the  2-D  analysis.  It  can  be  seen  that 
the  model  operates  as  expected,  with  the  injector  and  slotted  flap  providing  attached 
flow. 
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Figure  67.  3-D  Model  of  ADV3  Extended 


Figure  68.  3-D  Model  of  ADV4 


97 


99 


Figure  70.  Combustion  Efficiency  for  ADV4  Models 


;-D  Analysis  of  ADV5  Show  Attached  Flow 


The  three-dimensional  model  was  analyzed  next,  and  the  results  are  shown  in 
Figures  73  and  63.  The  mixing  performance  was  the  best,  by  far,  of  any  of  the  models 
tested,  reaching  100%  combustion  by  nearly  1/3  the  distance  of  the  ADV2  a20  model. 
The  pressure  loss  was  better  than  for  other  models  of  similar  blockage,  due  to  the 
net  gain  in  total  pressure  from  the  tangential  fuel  injection.  The  pressure  loss  is  still 
substantially  higher  than  the  baseline  IFF  pressure  loss  (by  a  factor  of  3),  and  the 
combustion  efficiency  gains  would  probably  not  overcome  the  pressure  losses. 

5.3.6  Advanced  Model  6 

Since  all  of  the  previously  tested  models  had  substantially  higher  pressure  losses 
than  the  baseline,  attention  turned  to  testing  models  with  smaller  base  areas  in  an 
effort  to  make  smaller  mixing  gains,  with  little  increase  in  pressure  loss.  The  first 
attempt  employed  a  small  triangular  fin  on  the  top  surface  of  the  IFF  near  the 
injector.  The  fin  measured  0.25  inches  high  and  0.4  inches  long,  and  was  placed  at  a 
30  degree  angle  to  the  flow.  The  leading  edge  was  placed  0.25  inches  outboard  of  the 
injector  centerline,  which  is  approximately  halfway  between  injectors.  ADV6  is 
shown  in  Figure  74,  with  the  crossflow  vectors  plotted  just  downstream  of  the 
flameholder.  The  analysis  showed  that  the  vortex  generator  produced  a  smaller, 
tighter  vortex  than  the  other  models,  but  has  a  pressure  comparable  to  the  baseline 
IFF  pressure  loss  (Figure  63).  The  numerical  test  was  isothermal  without  fuel 
injection,  to  examine  only  the  effects  of  the  fin.  To  judge  the  mixing  ability  of 
ADV6,  an  analysis  determined  that  a  fluid  element  passing  near  the  vortex 
generator  would  make  one  complete  turn  at  2.7  inches  downstream,  as  compared  to 
5  inches  downstream  for  ADV2  a20.  However,  the  vertical  extent  of  the  effect 
seemed  to  be  limited.  At  the  midpoint  between  the  IFF  surface  and  the  top  wall, 
there  was  very  little  influence  from  the  vortex  generator.  But  the  vortex  generator 
did  improve  the  lateral  fuel  distribution  along  the  base  of  the  flameholder. 

A  parametric  test  was  initiated  to  determine  the  effect  of  the  lateral  placement  of  the 
fin,  on  pressure  loss.  However,  when  the  fin  was  placed  near  either  symmetry 
plane,  numerical  instabilities  resulted,  possibly  from  unsteady  separation  from  the 
underside  of  the  fin. 
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Figure  73.  Combustion  Efficiency  for  ADV5  Models 


Figure  74.  3-D  Model  of  ADV6  with  Velocity  Vectors  Downstream  of  Flameholder 


5.3.7  Advanced  Model  7 

Another  approach  to  improving  the  performance  of  the  IFF,  was  to  take  advantage 
of  the  unsteady  vortices  that  were  shed  from  the  base  of  the  IFF  during  the 
isothermal  tests  (see  Section  4).  Advanced  model  ADV7  was  an  attempt  to  modify 
the  paths  of  the  shed  vortices.  Then,  by  altering  the  effect  laterally,  a  wave  could  be 
imposed  on  the  wake,  creating  a  larger  shear  surface  for  increased  mixing.  The  new 
design  had  a  simple  slant  cut  into  the  back  face  of  the  IFF  (Figure  75),  in  an  attempt 
to  skew  the  wake.  A  parametric  study  was  run  with  three  different  backface  angles. 
The  analysis  was  completed  in  two-dimensions  using  the  unsteady  solution 
algorithm.  While  each  geometry  produced  a  slightly  different  wake  pattern,  the 
vortices  simply  filled  any  voids  created  on  the  back  face  of  the  IFF,  with  no  change 
in  the  wake  overall.  Therefore,  the  design  was  not  successful.  Subsequent 
experimental  testing  verified  the  results. 


(a)  2-D  Model  of  ADV7  030 


(b)  2-D  Model  of  ADV7  a45 


(d)  3-D  Model  of  ADV7  a30 
Figure  75.  2-D  and  3-D  Models  of  ADV7 
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5.3.8  Advanced  Model  8 

The  final  advanced  model  tested  (ADV8)  again  used  the  momentum  of  the  fuel  to 
try  to  control  the  motion  of  the  wake.  Shown  in  Figure  76,  ADV8  followed  the 
baseline  IFF  shape,  but  injected  the  fuel  directly  into  the  wake  at  the  top  edge  of  the 
base  of  the  IFF.  Two  different  fuel  injection  angles  were  tested:  45  degrees,  and  67 
degrees  from  horizontal.  Again,  the  tests  were  conducted  in  two-dimensions  using 
the  unsteady  flow  solver. 


The  time-averaged  velocity  fields  are  shown  in  Figure  77  for  both  fuel  injection 
angles.  The  45  degree  case  did  exhibit  some  unsteady  vortex  motion  with  stronger 
perturbations  to  the  mean  flow  than  the  baseline  IFF.  However,  the  mean 
recirculation  zone  behind  the  IFF  was  considerably  smaller,  which  would  result  in 
short  residence  times  and  poor  flame  stability  characteristics.  The  67  degree  case  did 
not  show  any  signs  of  vortex  shedding,  and  therefore  offers  much  less  mixing  than 
the  baseline  IFF.  Since  the  fuel  is  injected  directly  into  the  wake,  the  mixing  of  the 
fuel  with  the  outer  regions  of  air  would  be  very  poor. 

5.4  Discussion 

Through  the  numerical  testing  of  the  many  different  flameholder  shapes, 
anchoring  the  flame  on  the  center  region,  as  opposed  to  the  ends  of  the  fingers, 
appears  preferable.  This  was  demonstrated  between  the  ADV2  a20  and  the  ADV3 
extended  model.  The  primary  difference  was  the  location  of  the  flameholder  base. 
ADV3  actually  had  more  geometric  blockage,  but  performed  better  in  terms  of 
combustion  efficiency  and  pressure  loss.  The  center  anchor  position  also  anchors 
the  flame  in  the  center  of  the  vortex.  At  lean  combustion  conditions,  with  the  onset 
of  vortex  shedding,  the  vortex  anchored  flame  may  prove  to  be  more  stable, 
reducing  the  lean  blowout  limit  of  the  combustor.  The  only  difficulty  with  the 
center  flameholder  comes  during  ignition.  Without  a  continuous  flameholder 
surface,  ignition  may  not  propagate  across  the  flameholder.  However,  the  fire  may 
travel  along  the  webs.  This  question  would  need  to  be  addressed  through 
experimental  testing. 
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The  best  design,  based  on  combustion  efficiency  and  pressure  loss,  was  the  ADV5 
model.  ADV5  provided  excellent  combustion  performance,  but  still  had  three  times 
the  pressure  loss  of  the  baseline.  Further  study  needs  to  be  completed  to  determine 
the  possible  tradeoffs  between  the  pressure  loss  and  the  mixing  performance. 
Possible  design  modifications  include  extending  the  Coanda  surface  and  eliminating 
the  flap.  The  vertical  extent  of  the  flameholding  surface  could  also  be  reduced.  The 
blockage  will  then  be  reduced,  but  the  fuel  will  still  be  providing  a  net  gain  to  the 
total  pressure  while  maintaining  a  strong  vortex.  Finally,  the  flameholding  surface 
could  be  reshaped  by  rounding  the  corners  and  allowing  it  to  follow  the  rounded 
shape  of  the  vortex  core. 

ADV5  has  another  problem  that  should  be  addressed.  The  design  relies  on  the  fuel 
to  provide  attached  flow  on  the  body.  The  combustor  would  not  perform  as  well 
when  operating  under  liquid  fuel  conditions,  or  under  low  power  conditions.  The 
extent  of  this  problem  has  not  been  determined. 

Design  ADV6  also  shows  some  promise,  with  the  use  of  simple  vortex  generators 
mounted  on  the  IFF  surface.  This  design  shows  the  best  promise  of  introducing 
some  performance  gain  with  very  little  additional  pressure  loss.  However,  the 
vortex  generators  would  be  in  immediate  danger  of  melting,  since  they  would  be 
immersed  in  an  auto-ignition  flame.  One  solution  may  be  to  inject  the  fuel  through 
the  trailing  edge  of  the  vortex  generators.  The  fuel  could  then  be  used  to  cool  the 
fins,  while  injecting  the  fuel  directly  into  the  vortex  core.  The  flameholder  and  the 
fins  would  then  be  out  of  the  fire,  with  a  corresponding  longer  operational  life.  In 
addition,  the  fuel  would  be  providing  a  gain  in  total  pressure  due  to  the  tangential 
injection.  This  can  be  a  detriment,  however,  as  the  fuel  may  not  penetrate 
svtfficiently  into  the  airstream.  The  vortex  generated  by  the  fins  is  not  as  large  as  the 
other  designs  (ADV2,  ADV5),  and  may  not  provide  enough  mixing  to  offset  the  loss 
of  fuel  penetration. 
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6.0  ISOTHERMAL  EXPERIMENTAL  TESTS:  ADVANCED  MODELS 

Four  advanced  models  were  fabricated  for  testing.  The  advanced  model,  ADV2, 
consisted  of  an  attachment  that  replaced  the  existing  0.125  inch  dump  plate  on  the 
baseline  EFF  (Figure  78).  The  attachments  consisted  of  three  identical  wedges  that 
created  45  degree  ramps  on  the  trailing  edge  of  the  IFF.  Each  wedge  was  0.5  inches 
long  and  extended  into  the  freestream  flow  by  0.5  inches.  The  separate  parts  were 
assembled  onto  a  0.125  inch  plate  that  was  then  fastened  to  the  baseline  IFF,  making 
the  total  length  of  the  ADV2  model  equal  to  1.5  inches.  Three  ramps  were 
fabricated,  with  the  center  ramp  facing  down,  and  the  two  outboard  ramps  facing  up. 
Only  three  ramps  were  included  in  the  tests  in  an  effort  to  simplify  the  numerical 
comparisons. 

Three  other  advanced  models  were  quickly  tested:  AD V7/ modi,  ADV7/mod2,  and 
ADV7/mod3  (see  Figure  79).  ADV7/modl  was  fabricated  by  taking  simple  trailing 
edge  extension  blocks  and  cutting  the  trailing  edge  back  at  a  30  degree  angle  to  the 
vertical.  The  model  consisted  of  two  extension  blocks,  one  cut  30  degrees  facing  up, 
the  other  cut  30  degrees  facing  down.  The  hope  was  that  each  piece  would  skew  the 
wake  in  a  favored  direction,  resulting  in  a  spanwise  pressure  gradient  that  would 
induce  a  vortex.  ADV7/mod2  attempted  to  further  enhance  the  mixing  by 
alternating  the  face  directions  across  the  span  of  the  model.  Finally,  ADV7/mod3 
took  the  mod2  extensions  and  mounted  them  on  the  leading  edge  of  the  baseline 
IFF,  replacing  the  hemispherical  nose. 

6.1  Test  Results 

The  only  advanced  model  extensively  tested  was  ADV2,  with  three  ramps.  The 
presence  of  an  axial  vortex  pair  was  evident  in  the  data,  with  a  peak  crossflow 
velocity  of  0.22  Uinf,  at  x/H  =  4.0.  A  parcel  of  air  crossing  near  this  point  would 
require  at  least  10  base  heights  to  make  one  full  revolution.  It  was  felt  by  the 
program  monitor  that  this  was  not  sufficient  to  substantially  enhance  the  fuel  air 
mixing. 

Another  model  tested,  ADV7/modl,  was  an  attempt  to  make  smaller  gains  at  fuel- 
air  mixing,  but  at  much  lower  pressure  losses.  However,  ADV7/modl  and  mod2 
both  produced  results  that  were  not  significantly  different  from  the  baseline  results, 
in  terms  of  both  the  velocities,  and  pressure  losses.  The  conclusion  was  that  the 
trailing  edge  can  be  modified  in  an  arbitrary  manner  with  no  effect  on  the  flowfield, 
as  long  as  no  protrusions  are  introduced  beyond  the  vertical  extent  of  the  baseline 
geometry. 
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a)  ADV7/mod  1  ;  two  spanwise  cuts 


b)  ADV7/mod  2  ;  multiple  spanwise  cuts 


c)  ADV7/mod  3  ;  two  spanwise  cuts  on  leading  edge 


Figure  79.  ADV7  Models  as  Tested 
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AD V7/ mods  did  produce  some  different  results.  The  goal  was  to  produce  leading 
edge  vortices  from  the  corners  of  the  leading  edge  extensions  as  the  flow  passed 
around  the  IFF.  Unfortunately,  no  corner  vortices  were  found,  although  a  large 
single  vortex  was  formed  through  the  entire  test  section  as  a  result  of  the  opposite 
direction  of  each  leading  edge  extension.  This  large  vortex  had  very  low  velocities 
and  would  not  significantly  enhance  mixing. 

6.2  Pressure  Losses 

The  static  pressure  losses  were  measured  for  all  models  based  on  the  difference 
between  the  static  pressure  at  x/H  =  -24,  and  x/H  =  48.  Both  pressure  ports  were 
connected  to  the  opposite  sides  of  a  micromanometer.  The  pressure  losses  were  low 
due  to  the  low  air  velocities  measuring  0.08%  for  the  baseline  IFF.  The  advanced 
model  ADV2  had  over  three  times  the  baseline  pressure  loss  at  0.29%.  This  value 
carmot  be  directly  compared  to  the  numerical  results,  since  the  fingers  did  no  extend 
to  the  walls.  Based  on  these  results,  it  was  felt  that  the  ADV2  pressure  losses  would 
offset  any  gains  made  from  enhanced  fuel /air  mixing  created  by  the  vortex 
formation.  All  of  the  ADV7  models  had  the  same  pressure  loss  as  the  baseline, 
since  they  all  had  the  same  base  areas  and  geometric  blockages  as  the  baseline. 
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7.0  EXPERIMENTAL  COMBUSTION  TESTS 


The  original  Phase  I  plans  called  for  the  design  and  fabrication  of  a  water-cooled 
combustion  test  section  for  installation  into  Test  Cell  19.  The  tunnel  design  placed 
the  IFF  model  spanning  the  tunnel  in  the  vertical  direction,  providing  a  top  view  of 
the  IFF  through  the  side  window.  The  overall  design  included  full  view  windows 
on  the  top  and  bottom  of  the  test  section,  and  partial  view  windows  on  the  sides. 
The  test  section  was  42  inches  long  and  the  optical  quality  quartz  windows  were  14 
inches  long.  The  window  retainers  could  be  reversed  to  change  the  window 
viewing  area  from  upstream  to  downstream,  providing  a  total  viewing  area  of  26 
inches  in  the  axial  direction. 

LDV  velocity  data  and  CARS  temperature  data  were  to  be  taken  throughout  the  test 
section  for  the  baseline  and  one  advanced  model,  at  one  equivalence  ratio.  The 
isothermal  tests  were  also  to  be  repeated  on  the  baseline  model.  Test  Cell  19 
included  support  for  providing  clean  nonvitiated  air  up  to  temperatures  of  1000  F, 
and  had  a  full  settling  chamber  and  transition  section  for  high  quality  inlet  flow.  In 
addition  to  a  new  test  section,  a  new  Mach  control  section  was  needed  for  the 
subsonic  tests. 

The  design  for  both  the  test  section  and  Mach  control  section  was  completed  on 
schedule.  Most  of  the  mechanical  design  and  drafting  work  was  completed  by  Mr. 
Gary  Haines  of  Cincinnati  Control  Dynamics,  Inc.  Mr.  Steven  Green  of  Wright 
Laboratory  completed  many  of  the  heat  transfer  and  cooling  requirement 
calculations  for  the  test  section.  However,  the  sections  were  never  installed  and  the 
tests  were  canceled  due  to  scheduling  difficulties. 

For  documentation  purposes,  the  basic  component  design  criteria  are  presented  in 
the  next  sections  for  each  component,  followed  by  a  list  of  uncompleted  items  for 
future  reference. 

7.1  Test  Section  Design 

The  conceptual  design  of  the  test  section  (see  Figure  80)  was  48  inches  long  with  a 
6x6  inch  internal  cross  section.  The  reversible  window  panels  provided  access  from 
6  inches  upstream  of  the  model  trailing  edge  to  20  inches  downstream,  on  three 
sides.  The  bottom  window  was  limited  to  view  from  the  trailing  edge  to  8.5  inches 
downstream.  All  solid  walls  of  the  tunnel  were  constructed  from  a  copper/ stainless 
steel  sandwich  with  cooling  channels  cut  into  the  copper.  The  outer  stainless  jacket 
sealed  against  the  copper  to  close  the  water  channels.  Water  channels  were  also 
provided  along  the  side  window  frames  to  provide  some  cooling  to  the  copper 
frame,  to  reduce  the  linear  thermal  growth  of  the  copper  wall. 
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Combustion  Test  Facility 

Cooling  Water  System 
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Figure  80.  Schematic  of  Combustion  Test  Section 


The  quartz  windows,  being  uncooled,  were  the  limiting  structural  components.  The 
maximum  overall  equivalence  ratio  was  limited  to  <j)=0.2  to  keep  the  window 
temperatures  below  1300  K.  Window  growth  with  temperature  was  not  considered 
a  problem  (<0.005  inch)  and  was  neglected.  A  window  thickness  of  0.75  inches  was 
chosen  for  all  windows.  Table  3  shows  the  window  dimensions  and  factor  of  safety 
for  all  windows  at  a  2  atmosphere  pressure  differential. 

Table  3.  Factor  of  Safety  for  Timnel  Windows  at  30  psig 


window  size  (inches) 

top  14x6.75 

sides  14x3 

bottom  8.5x6.75 


factor  of  safety 
4 
17 
7 


The  side  windows  were  held  in  place  by  retaining  strips.  The  side  windows  become 
part  of  the  side  wall  window  retainer  unit  when  installed.  The  top  and  bottom 
windows  do  not  become  part  of  their  respective  wall  units,  but  the  windows  are 
trapped  between  the  stainless  steel  outer  jacket  and  the  edges  of  the  copper  side 
walls.  Therefore,  if  the  top  wall  window  retainer  is  removed,  the  window  will 
remain  resting  on  top  of  the  test  section,  supported  by  the  side  walls. 

To  facilitate  the  acquisition  of  tunnel  temperatures  and  pressures,  taps  were  drilled 
in  several  places  down  the  midline  of  each  wall  section.  The  ports  open  up  to  a 
larger  threaded  chamber  in  the  stainless  steel,  and  should  be  sealed  with  a  pipe  plug 
when  not  in  use. 


The  test  section  was  fabricated  at  Sydney  Tool  and  Die  in  Sydney,  Ohio.  The 
components  were  shipped  to  Wright  Patterson  AFB,  and  inspected.  Each  wall  was 
assembled  and  pressure  checked  to  100  psi.  The  windows  were  cut  and  polished  by 
Glass  Fab,  Inc,  located  in  Rochester,  N.Y.  The  specifications  on  the  windows 
required  an  80/50  scratch/dig  polish,  a  flatness  of  1  fringe  per  inch,  and  a 
homogeneity  of  5x1 0-^  or  better. 

7.2  Transition  Section  Design 

When  the  tests  were  rescheduled  for  Test  Cell  18,  a  new  transition  had  to  be 
designed  to  fit  the  new  facility.  The  transition  had  one  main  purpose  -  to  smoothly 
transition  from  the  round  pipe  of  the  tunnel  air  supply,  to  the  square  test  section. 
To  reduce  costs,  standard  hardware  was  utilized  as  much  as  possible. 

The  critical  component  of  the  entire  section  was  the  transition  liner,  and  was 
custom  made.  The  liner  internal  geometry  was  designed  using  a  purely 
mathematical  function  to  transition  from  a  round  cross  section  to  a  square  cross 
section  with  zero  axial  slope  at  each  end.  36  The  transition  function  had  no  physical 
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basis.  The  liner  shape  was  then  electronically  delivered  to  the  Air  Force  Machine 
Shop  at  Wright  Laboratory,  where  stereo  lithography  was  used  to  fabricate  the  liner. 
Typically  used  to  fabricate  manufacturing  prototypes,  stereo  lithography  uses  an 
ultra  violet  laser  to  paint  a  given  cross  section  of  the  design  in  a  thin  layer  of  liquid 
polymer.  The  liquid  polymer  hardens  when  struck  with  the  laser.  The  liquid 
polymer  level  is  then  raised  slightly  in  the  tank,  and  the  next  cross  section  is 
painted.  The  process  is  repeated,  building  up  one  cross  section  on  top  of  the  other 
until  the  article  is  completed.  There  is  no  constraint  to  the  complexity  of  the  article 
that  can  be  produced  with  stereo  lithography. 

The  transition  liner  was  actually  fabricated  in  three  sections,  to  make  the  most 
efficient  use  of  the  stereo  lithography  machine.  When  finished,  the  pieces  were 
sanded  and  glued  together  to  form  the  finished  product.  Although  the  process  is 
typically  used  to  construct  prototypes  or  mold  casts,  the  plastic  transition  liner  was 
the  actual  part  used  in  the  ttmnel.  The  liner  acted  only  to  provide  a  flow  path,  and 
did  not  carry  any  loads  or  contain  any  pressure.  To  contain  the  pressure,  the  liner 
was  fit  inside  of  a  12-inch  nominal  pipe,  that  was  welded  closed  at  one  end,  and 
fitted  with  an  open  flange  at  the  other.  The  closed  end  was  machined  with  a  rotmd 
hole  to  hold  the  round  end  of  the  liner,  and  machined  with  threaded  holes  to  mate 
to  the  existing  tunnel  supply  pipe.  The  flange  end  was  fit  with  a  blank  flange  that 
had  been  machined  with  a  square  hole  to  fit  the  square  end  of  the  liner.  The  blank 
flange  also  had  threaded  holes  to  connect  to  the  test  section.  To  protect  the  liner 
from  vibrational  stress,  and  to  provide  additional  support,  expanding  foam  filler 
was  purchased  to  fill  the  space  between  the  pipe  and  the  outside  of  the  liner. 
However,  this  step  was  never  completed. 

7.3  Choke  Block  Design 

The  design  for  the  choke  block  section  did  not  change  when  the  tests  were  moved 
from  Test  Cell  19  to  Test  Cell  18.  Several  designs  were  considered,  including  a  water 
cooled  cone  that  moved  axially  to  change  the  blockage  at  the  test  section  exit.  While 
the  cone  provides  the  easiest  method  of  tunnel  operation,  simple  2-D  flaps  were 
selected  due  to  the  ease  of  design  and  fabrication. 

The  choke  blocks  were  designed  with  a  different  approach  than  the  test  section. 
Instead  of  a  copper  cooling  wall  covered  by  a  stainless  steel  jacket,  the  choke  blocks 
used  single  plate  stainless  steel  walls  with  straight  axial  cooling  tubes  drilled 
through  the  ends  of  the  top  and  bottom  walls.  The  water  flowed  from  the  front  of 
the  section  and  dumped  out  the  rear  of  the  walls  directly  into  the  tunnel  exhaust 
section.  The  side  walls  were  not  cooled  since  they  were  protected  by  the  flappers. 
The  copper  flappers  were  cooled  by  water  sprayed  on  the  leading  edge  from  holes  in 
the  side  walls,  just  upstream  of  the  flappers.  The  Mach  number  was  controlled  by 
setting  the  position  of  each  flapper  with  a  threaded  plunger  mounted  near  the 
trailing  edge  of  each  flapper.  Care  must  be  taken  to  set  each  flapper  symmetrically. 
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or  flow  asymmetries  could  develop  in  the  test  section.  The  section  was  fabricated  by 
Sydney  Tool  and  Die. 

7.4  Baseline  Model  Design 

The  baseline  model  designed  for  the  combustion  tests  was  less  complex  than  the 
isothermal  model  design.  For  the  isothermal  model,  several  tail  pieces  were 
designed  to  fit  on  the  centerbody  in  a  modular  fashion,  and  without  any  screw  heads 
exposed  on  the  back  face.  The  combustion  model  design  was  not  as  modular.  The 
model  was  machined  from  a  solid  piece  of  stainless  steel  to  the  shape  of  the  IFF, 
with  an  additional  tail  piece  that  screwed  directly  onto  the  base  of  the  primary  piece. 
The  screw  heads  were  countersunk  to  provide  minimum  flow  interference.  The 
TrO-ittchr^ase  height  of  the  model  provided  the  same  blockage  as  in  the  isothermal 
tests.  A  pair  of  holes  were  drilled  through  the  span  of  the  model  to  serve  as  fuel 
manifolds.  Stainless  steel,  thick-walled  tubes  were  threaded  into  each  hole,  passing 
through  the  floor  of  the  tunnel,  and  acted  as  two  pins  to  maintain  the  model's 
orientation  in  the  tunnel.  To  hold  the  model  firmly  against  the  tunnel  floor,  two 
snap  rings  were  installed  near  the  midpoint  of  the  length  of  the  fuel  tubes,  capturing 
a  small  steel  plate  between  the  rings  and  the  outside  of  the  tunnel  floor  (Figure  81). 
Four  bolts  threaded  through  the  plate  and  pressed  against  the  outside  floor  of  the 
tunnel.  By  uniformly  tightening  the  bolts,  the  plate  exerted  a  uniform  pressure 
against  the  snap  rings,  firmly  holding  the  IFF  model  in  place.  The  opposite  end  of 
the  model  was  fit  with  an  o-ring  to  provide  a  cushion  against  the  upper  window. 
The  Model  was  fabricated  by  Tipp  Machine  &  Tool,  Inc,  in  Tipp  City,  Ohio. 

7.5  Miscellaneous  Information 

While  many  of  the  parts  have  been  assembled,  the  entire  system  has  never  been 
assembled,  and  no  leak  checks  have  been  completed  on  the  air  circuit.  Several  sub¬ 
components  need  to  be  obtained  or  fabricated  to  complete  the  system.  These 
include: 

1)  leveling  bolts; 

2)  gaskets  between  components; 

3)  choke  block  set  screws; 

4)  retaining  snap  rings  for  IFF  model;  and 

5)  installation  of  foam  filler  in  transition  section. 

For  the  leveling  bolts,  two  1/2  inch  bolts  should  be  sufficient  on  each  end  of  the  test 
section.  The  bolts  will  be  threaded  through  the  floor  of  the  mounting  structure  and 
rest  against  the  bottom  of  the  tunnel. 
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Figure  81.  Schematic  of  IFF  Mounting  System  for  Combustion  Tests 


Gaskets  will  need  to  be  fit  between  the  following  components: 


a)  air  supply  pipe; 

b)  transition  section; 

c)  mounting  flange; 

d)  test  section; 

e)  choke  block  section; 

f)  downstream  mounting  flange;  and 

g)  tunnel  exhaust. 


The  choke  block  set  screws  will  have  to  be  made  from  1/2  inch  fine  thread  all¬ 
thread,  and  should  be  at  least  4  inches  long.  A  nut  tack-welded  on  the  end  of  the  all¬ 
thread  should  provide  as  a  suitable  bolt  head.  During  the  design  process,  it  was 
suggested  that  the  set  screws  be  fabricated  with  a  micrometer-type  scale  to  provide 
precise  settings  on  both  sides  of  the  tunnel.  It  was  decided  that  such  accuracy  was 
not  necessary.  A  built-in  device  to  measure  the  distance  from  the  side  of  the  tunnel 
with  an  accuracy  of  0.050  inch  should  be  sufficient. 


The  snap  rings  used  to  hold  the  IFF  model  in  place  still  need  to  be  acquired.  The 
snap  rings  are  an  odd  size  and  may  have  to  be  special  ordered  (Truarc  external  series 
5100-40, 13/38) 


The  hmnel  design  was  carefully  considered  and  well  constructed.  However,  there 
were  some  lessons  learned  during  the  fabrication  process  that  should  be  considered 
for  future  modifications.  The  top  window  retainer  should  be  held  in  by  positive 
locking  latches  to  make  removal  of  the  window  easier.  Only  one  window  needs 
easy  removal,  to  provide  access  to  the  model  and  the  other  windows  for  cleaning. 
In  addition,  all  window  retainers  need  handles  installed  for  easier  handling.  The 
retainers  fit  with  tight  tolerances  and  are  difficult  to  remove.  Without  handles,  the 
chances  are  increased  that  a  section  will  be  dropped  during  installation  or  removal. 

7.6  Test  Cancellation 

The  combustion  tests  were  originally  scheduled  to  be  completed  in  Test  Cell  19,  and 
were  to  use  the  clean  heated  air  system  available  in  that  facility.  Scheduling 
overruns  of  other  tests  operating  in  Cell  19  at  the  time  forced  the  current  project  to 
be  rescheduled  for  Test  Cell  18.  The  rescheduling  had  two  major  implications: 
1)  there  was  no  heated  air  system  available  in  Cell  18,  and  2)  the  new  test  section 
would  have  to  be  retrofitted  to  fit  the  new  facility,  including  the  fabrication  of  the 
transition  section.  Both  were  deemed  acceptable  to  keep  the  tests  on-line. 
Unfortunately,  there  were  also  other  tests  being  conducted  in  Test  Cell  18.  The 
current  tests  were  postponed  a  number  of  times,  and  were  eventually  canceled 
when  the  program  ran  out  of  time. 
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8.0  NUMERICAL  COMBUSTION  TESTS 


Since  no  experimental  combustion  data  were  produced  in  this  program,  another 
data  set  was  selected  for  the  numerical  study.  The  data  set  selected  for  study  was  that 
of  Sjunnesson,  et  al.37  This  data  set  was  selected  because  it  was  a  recent  experiment 
that  included  detailed  CARS,  LDV,  and  gas  analyses  for  a  simple  geometry.  The  data 
set  had  been  the  focus  of  several  earlier  numerical  studies,  and  was  judged  to  be  a 
good  data  base  for  CFD  comparisons. 

8.1  Experimental  Data  Set 

The  flameholder  model  tested  by  Sjunnesson  was  cylindrical  with  an  equilateral 
triangular  cross  section  measuring  40  mm  on  a  side  (see  Figure  82).  The  tunnel 
cross  section  was  12  cm  high  by  24  cm  wide.  A  premixed  air/ propane  mixture  was 
introduced  to  the  tunnel  at  equivalence  ratios  varying  from  0  to  1.1.  Inlet  Mach 
numbers  ranged  from  0  to  0.2,  at  a  Reynolds  number  of  5.5x105,  based  on  the 
hydraulic  diameter  of  the  combustor.  Inlet  temperature  was  288  K,  and  the  pressure 
was  ambient. 


Figure  82.  A  Schematic  Drawing  of  the  Validation  Rig 
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The  isothermal  tests  showed  expected  results.  Vortices  were  shed  behind  the 
flameholder  in  a  regular  pattern.  At  the  combusting  conditions,  no  vortex  shedding 
was  indicated,  even  at  the  lowest  equivalence  ratio  (d=0.58).  Tunnel  acoustics  are 
known  to  play  an  important  role  in  ducted  burners,  and  the  experimentalists 
noticed  three  instabilities  during  their  tests.  A  screech  mode  (1400  Hz)  was  found  at 
(j)=0.72,  a  longitudinal  mode  (100  Hz)  was  found  at  (|)  =  0.95,  and  a  combination  of  the 
two  modes  was  found  at  the  same  conditions  as  the  screech  mode.  The  low 
frequency  mode  can  cause  numerical  difficulties,  since  it  is  of  the  same  frequency 
range  as  the  vortex  shedding,  and  can  have  significant  amplitudes.  Sjunnesson 
found  that  the  low  frequency  mode  created  large  changes  in  the  structure  of  the 
wake  near  the  flameholder.  At  one  point  in  the  periodic  cycle,  the  wake  bulged  out 
and  became  much  thicker  than  the  base  height  of  the  flameholder. 

Numerical  Details 

The  CFD  code  tested  with  the  triangular  flameholder  was  the  same  code  used  for  the 
isothermal  tests:  CFD-ACE.  CFD- ACE  has  several  combustion  models,  including 
simple  instantaneous  chemistry,  equilibrium  chemistry,  and  finite-rate  chemistry 
using  one,  two,  or  four  step  models.  The  finite-rate  methods  can  be  made  to  react  to 
equilibrium  conditions  using  an  efficient  table  look-up  routine.  Several  advanced 
methods  are  also  available  in  CFD-ACE,  such  as  the  prescribed  PDF  and  Monte  Carlo 
PDF  models. 

The  numerical  tests  were  completed  using  the  single  step  propane  model  with  the 
reaction  rate  constants  taken  from  Westbrook  and  Dryer^^; 

A  =  8.6  X  ion  a  =  0.1 

n  =  0  b  =  1.65 

Ea  =  30.0 

The  simulations  were  performed  at  the  following  freestream  conditions: 

Uco  =  17.3  misec 
=  100  kPa 

T==  =  288K 

(p  =  0.61 


The  grid  is  shown  in  Figure  83.  Two  different  grids  were  used,  with  the  difference 
coming  in  the  location  of  the  inlet  plane.  The  first  grid  started  at  200  mm  upstream. 
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and  used  the  first  plane  of  experimental  data  for  the  inflow  conditions.  However, 
due  to  interactions  with  the  tunnel  acoustics,  it  was  felt  that  the  entire  tunnel  length 
had  to  be  modeled,  and  the  grid  was  started  at  819  mm  upstream  of  the  flameholder, 
with  a  uniform  inlet  condition.  No-slip  walls  were  maintained  on  both  walls  of  the 
tunnel  and  on  the  walls  of  the  flameholder.  The  downstream  grid  exited  into  a 
dump  diffuser  with  an  equivalent  2-D  area  as  the  experimental  axisymmetric  dump 
chamber.  The  exit  plane  of  the  dump  chamber  was  held  at  a  constant  pressure 
condition. 

The  numerical  tests  were  completed  using  the  same  methodology  as  used  for  the 
isothermal  unsteady  analyses.  Therefore,  any  vortex  shedding,  or  lack  of  vortex 
shedding,  along  with  other  instabilities,  would  be  accounted  for  in  the  solution. 

8.3  Numerical  Results 

The  case  selected  for  study  was  at  (p  =  0.61,  corresponding  to  case  1  from  Sjimnesson 
et  al.37  The  objective  of  the  calculation  was  to  compare  numerical  results  with 
experimental  measurements.  The  experimental  study  did  not  show  any 
combustion  instabilities  at  this  equivalence  ratio.  High  speed  photographs  indicated 
the  heat  release  occiirred  in  the  simultaneous  roll-up  of  the  two  shear  layers  behind 
the  flameholder.  The  high  strain  rate  reduced  the  heat  release  immediately  behind 
the  flameholder.  Temperature  measurements  on  the  combustor  centerline  showed 
unburned  fluid  was  entrained  behind  the  recirculation  zone. 

The  numerical  results  showed  strong  pressure  oscillations,  which  were  not  expected 
at  this  equivalence  ratio.  These  pressure  oscillations  appeared  to  be  an  axial  "organ- 
pipe"  mode  with  a  frequency  of  approximately  200  Hz.  This  mode  is  similar  to  the 
combustion  instability  observed  experimentally  at  a  higher  equivalence  ratio.  The 
pressure  history  at  a  point  on  the  combustor  centerline,  behind  the  flameholder,  is 
shown  in  Figure  84.  The  pressure  throughout  the  combustor  oscillated  in  phase, 
with  the  maximum  amplitude  at  the  inlet. 

The  flame  shape  in  the  region  immediately  behind  the  flameholder  is  shown  at  a 
sequence  of  different  times  in  Figure  85.  The  difference  in  time  between  plots  is 
0.4  ms.  The  pressure  at  the  base  of  the  flameholder  is  at  a  maximum  at  time  tj. 
Shortly  after  the  pressure  reaches  the  peak,  a  pair  of  counter-rotating  vortices  are 
shed  simultaneously  from  the  corners  of  the  flameholder  and  the  flame  front  bulges 
out.  The  pressure  peaks  again  at  fs,  triggering  the  shedding  of  another  vortex  pair. 
As  the  vortices  are  convected  downstream,  cold  fluid  is  entrained  into  the  wake 
region.  The  flow  becomes  asymmetrical  further  downstream,  as  the  vortices 
interact.  This  flow  did  not  reach  a  regular  periodic  state. 
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Figure  83.  Grid  for  Reacting-Flow  Case  Modeling  Sjunnesson  Experiment 
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Figure  84.  Pressure  History  at  Centerline,  x/H  =  1,  Demonstrates  Acoustic  Pressure 
Waves  in  Tunnel 


Figure  85.  Sequence  of  Flame  Shapes  for  Triangular  Flameholder  with  (p  =  0.61; 
Color  Denotes  Temperature  (Blue  =  288  K,  Red  =  1710  K) 
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Effects  of  Initial  and  Boundary  Conditions 


The  pressure  oscillations  that  dominated  the  flow  were  created  by  the  initial 
conditions  used  for  the  transient  combustion  simulation,  and  then  sustained  by 
interactions  with  the  heat  release  from  combustion  during  the  simulation.  The 
initial  conditions  for  the  flow  were  uniform  and  the  same  as  the  upstream 
boundary,  except  for  the  region  behind  the  flameholder.  A  high  temperature  was 
specified  behind  the  flameholder  to  start  the  finite-rate  reaction.  This  is  the  same 
methodology  used  for  the  nonreacting  calculations,  except  for  the  elevated 
temperature.  These  conditions  are  artificial  and  do  not  satisfy  conservation  of  mass, 
momentum,  or  energy.  After  a  transition  period,  the  initial  conditions  are  not 
expected  to  be  important.  The  flow  rate  at  the  inlet  and  the  pressure  at  the  exit  are 
fixed,  so  that  unsteady  flow  patterns  should  be  caused  only  by  instabilities  in  the 
flow. 

Integration  of  the  continuity  equation  shows  the  local  rate  of  change  of  density  is 
proportional  to  the  mass  imbalance  for  each  control  volume.  The  continuity 
equation  is  used  to  derive  a  pressure  correction  equation  for  the  numerical  method 
used  in  this  study.  The  mass  imbalance  for  each  cell  is  part  of  the  source  term  for 
the  pressure  correction  equation.  The  mass  imbalances  caused  by  the  initial  flow 
field  not  satisfying  continuity  create  large  initial  changes  in  pressure.  These  changes 
were  expected  to  decay  to  a  much  lower  level. 

The  boxmdary  conditions  for  this  simulation  reflect  some  of  the  pressure  wave  that 
is  generated  by  the  initial  conditions.  The  pressure  at  the  exit  boundary  was  fixed, 
while  the  pressure  at  the  inlet  boundary  had  a  zero  gradient.  As  a  result,  a  standing 
quarter-wave  was  set  up  in  the  combustor.  Changes  in  pressure  led,  in  turn,  to 
changes  in  enthalpy,  which  affected  the  local  reaction  rate  and  changed  the  flame 
shape.  The  flame  shape  was  also  altered  by  changes  in  the  local  flow  patterns  caused 
by  pressure  oscillations.  Unsteady  heat  release  causes  combustion  instability  when 
the  fluctuating  pressure  is  in  phase  with  the  fluctuating  pressure  field.  The 
unsteady  heat  release  appeared  to  resonate  with  the  pressure  field  in  this  simulation 
and  sustain  the  pressure  oscillations. 

One  method  used  to  reduce  the  transition  period  at  the  beginning  of  the  transient 
simulation  and  the  disturbance  in  the  pressure  field  was  to  first  perform  a  steady- 
state  calculation.  The  steady  state  results  would  be  used  as  the  initial  conditions  for 
the  transient  case.  This  procedure  was  used  in  this  study.  The  residuals  for  each 
variable  could  only  be  reduced  two  to  three  orders  of  magnitude  before  the  solution 
began  to  oscillate.  The  solution  was  not  symmetric  and  had  the  appearance  of  the 
alternate  vortex  shedding  seen  in  the  nonreacting  simulations.  The  partially- 
converged  steady-state  solution  was  then  used  as  the  initial  value  for  the  transient 
calculation.  Errors  in  mass  conservation  were  reduced,  but  not  eliminated,  by  the 
steady-state  calculation.  The  mass  conservation  errors  were  still  large  enough  to 
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create  large  pressure  oscillations.  There  was  no  substantial  difference  in  the  flow 
patterns  in  simulations  started  from  steady  state  results  or  from  the  initial 
conditions  discussed  above. 

Flat  Plate  Flameholder 

The  flow  for  this  case  was  expected  to  show  vortex  shedding  patterns  similar  to  the 
nonreacting  cases,  although  the  volumetric  expansion  caused  by  heat  release  from 
the  reaction  reduces  the  vortidty  of  the  fluid  and  is  expected  to  stabilize  the  shear 
layers.  An  earlier  study  of  unsteady  combustions^  showed  the  transition  process  as 
the  equivalence  ratio  was  reduced  for  a  flame  stabilized  by  a  flat  plate.  The  operating 
conditions  for  the  earlier  study  used  a  higher  inlet  temperature  (800K)  and  lower 
pressure  (30  kPa),  but  the  same  fuel  and  reaction  mechanism.  The  flow  field  was 
steady  as  the  equivalence  ratio  was  reduced  from  1.0  to  0.8.  As  the  equivalence  ratio 
was  further  reduced  to  0.4,  low-amplitude,  high-frequency  oscillations  were 
observed  in  the  pressure  field  but  the  flow  field  remained  steady.  As  the 
equivalence  ratio  was  reduced  even  further,  vortices  were  shed  from  the  top  and 
bottom  of  the  flameholder  in  an  alternating  fashion.  Low-frequency  pressure 
oscillations  were  observed  upon  each  change  in  the  inlet  fuel /air  ratio.  These 
oscillations  died  out  after  several  cycles.  A  regular,  periodic  flow  was  observed  for 
the  lower  equivalence  ratios. 

There  was  some  concern  that  the  difference  in  flow  patterns  observed  in  this  study 
and  the  earlier  study  with  the  flat  plate  were  due  to  changes  in  the  solution 
algorithms  during  the  development  of  the  computer  code.  Significant  changes  had 
been  made  to  the  code,  including  the  addition  of  multiblock  capability  and  major 
changes  to  the  solution  of  the  energy  equation.  In  an  effort  to  determine  if  code 
changes  were  responsible  for  the  unexpected  flow  patterns,  one  of  the  flow 
conditions  that  produced  vortex  shedding  in  the  earlier  study  was  reproduced  The 
previous  results  were  successfully  reproduced  with  the  computer  code  used  in  this 
study,  indicating  the  differences  in  the  observed  flow  patterns  were  due  to  physical 
differences  in  the  flows. 

One  significant  difference  between  this  and  the  earlier  study,  besides  the  operating 
conditions,  was  that  the  earlier  study  used  a  constant  turbulent  viscosity,  instead  of 
the  two-equation  model  used  here.  In  order  to  determine  that  the  different  flow 
pattern  was  caused  by  physics,  and  not  modeling,  a  calculation  was  made  with  a  flat 
plate  flameholder,  constant  turbulent  viscosity,  and  the  operating  conditions  for  the 
triangular  flameholder.  The  plate  height  was  the  same  as  the  flameholder  base  and 
the  tunnel  dimensions  were  the  same  as  used  for  the  triangular  flameholder. 
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Results  with  the  flat  plate  were  consistent  with  the  triangular  flameholder  case.  The 
flow  was  nearly  symmetrical  and  dominated  by  axial  pressure  oscillations. 
Instantaneous  velocity  vectors  and  temperature  contours  are  shown  for  the  region 
immediately  behind  the  flameholder  in  Figure  86.  The  pair  of  vortices  are 
periodically  shed  and  convected  downstream.  Unburnt  gas  is  entrained  by  the 
vortex  shedding  and  intermittently  reaches  the  centerline.  Pressure  and 
temperature  histories  on  the  centerline  are  shown  in  Figure  87.  The  temperature 
drop  when  unreacted  fluid  reaches  the  centerline  is  strongly  correlated  with  peaks 
in  the  pressure  field,  although  the  cold  fluid  does  not  reach  the  centerline  on  every 
cycle.  The  similarity  between  the  results  with  both  flameholder  shape  shows  the 
flow  patterns  are  not  an  artifact  of  the  turbulence  modeling  or  multiblock  solution. 

Because  the  experimental  instabilities  were  only  observed  for  higher  equivalence 
ratios,  the  calculation  was  repeated  with  the  equivalence  ratio  reduced  to  0.55. 
Instantaneous  velocity  vectors  and  temperature  contours  are  shown  in  Figure  88  for 
the  same  region  shown  above.  The  vortices  are  now  shed  alternately  from  the  top 
and  bottom  of  the  flameholder.  Pressure  oscillations  are  an  order  of  magnitude 
smaller  than  with  9  =  0.61.  The  reduced  amount  of  heat  released  at  this  equivalence 
ratio  is  not  sufficient  to  sustain  the  pressure  oscillations,  and  the  flow  is  controlled 
by  the  instability  of  the  shear  layers  behind  the  flameholder. 

Reduced  Equivalence  Ratio 

Because  the  flat  plate  results  showed  reduced  pressure  oscillations  at  a  lower 
equivalence  ratio,  the  triangular  flameholder  calculation  was  repeated  with  an 
equivalence  ratio  of  0.55.  All  conditions  except  the  equivalence  ratio  were 
imchanged.  This  calculation  did  not  show  the  alternate  vortex  shedding  observed 
for  the  flat  plate  flameholder  with  the  same  conditions.  Figure  89  shows  velocity 
vectors  and  temperature  contours  in  the  vicinity  of  the  flameholder.  The  flow 
pattern  is  symmetric.  There  is  a  small  oscillation  of  the  wake  region,  but  no 
evidence  of  large-scale  vortex  shedding. 

The  simulation  was  started  with  the  RNG  k-e  turbulence  model.  The  final  5000 
time  steps  were  made  with  a  constant  effective  viscosity,  to  check  the  sensitivity  of 
the  solution  to  the  turbulence  modeling.  The  history  of  pressure  behind  the 
flameholder  is  shown  in  Figure  90,  with  time  t  =  0.0  corresponding  to  the  change  in 
turbulence  model.  The  amplitude  of  the  pressure  oscillations  decays,  indicating  the 
reduced  heat  release  does  not  sustain  the  axial  instability  mode.  The  amount  of  heat 
release  appears  to  be  enough  to  prevent  the  alternating  vortex  shedding.  The 
difference  in  flow  patterns  between  the  flat  plate  and  triangular  flameholders  must 
be  attributed  to  the  change  in  the  shape  of  the  flameholder. 
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Figure  86.  Instantaneous  Velocity  Vectors  and  Temperature  Contours  for  the  Flat 
Plate  Flameholder  with  cp  =  0.61 


(a)  (b) 


Figure  87.  Time  Histories  of  Pressure  and  Temperature  Behind  the  Flat  Plate 
Flameholder  (cp  =  0.61) 


Figure  88.  Instantaneous  Velocity  Vectors  and  Tempers  iure  Contours  for  Flat 
Plate  Flameholder  with  cp  =  0.55 


Figure  89.  Instantaneous  Velocity  Vectors  and  Temperature  Contours  for 
Triangular  Flameholder  with  (p  =  0.55 
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t 

Figiare  90.  Pressure  Piistory  Behind  Triangular  Flanneholder  with  (p  =  0.55 
Summary 

No  comparison  has  been  made  to  experimental  results  because  the  flow  regime  does 
not  appear  to  be  the  same  at  the  equivalence  ratio  considered.  The  occurrence  of  an 
axial  combustion  instability  mode  at  a  lower  equivalence  ratio  than  was  seen 
experimentally  could  be  due  to  several  factors.  First,  the  boimdary  conditions  in  the 
simulation  may  be  artificially  reflecting  pressure  waves  back  into  the  computational 
domain.  Care  was  taken  to  locate  the  boundaries  at  the  combustor  inlet  and  exit 
planes,  but  information  is  transferred  into  the  computational  domain  at  the  exit  by  a 
characteristic  wave  moving  upstream. 

A  second  factor  in  the  prediction  of  combustion  instability  is  the  combustion  model. 
A  single-step  model  was  used  in  this  study.  While  this  model  was  successful  in 
earlier  studies  and  with  simpler  geometry,  the  immediate  conversion  of  fuel  into 
products  could  alter  the  way  heat  release  from  the  reaction  interacts  with  the 
fluctuating  pressure  field.  A  more  realistic  model  would  have  multiple  steps  and 
the  heat  wo^d  not  be  released  all  at  once. 

A  third  factor  is  the  interaction  between  turbulence  and  combustion.  The 
combustion  model  used  in  this  study  is  "quasi-laminar"',  in  that  the  production  rate 
is  calculated  from  the  mean  species  concentrations.  The  experimental  results 
showed  the  heat  release  was  controlled  by  mixing  in  the  shear  layers.  In  the 
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numerical  study,  the  reaction  rate  was  fast  enough  for  the  fuel  to  completely  react  at 
the  base  of  the  flameholder. 

More  work  needs  to  be  done  on  the  transient  combustion  simulations  to  distinguish 
between  numerical  and  physical  instabilities.  The  initial  and  boundary  conditions 
must  be  as  physically  correct  as  possible.  The  combustion  simulations  should  be 
performed  with  more  realistic  chemistry  models  that  account  for  the  effect  of 
turbulence  on  combustion. 
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9.0  CONCLUSIONS 


Numerical  analyses  using  2-D  and  3-D  steady-state  and  time  accurate  methods, 
under  reacting  and  nonreacting  conditions,  were  performed  on  the  baseline  IFF 
geometry  and  28  advanced  designs.  Many  of  the  advanced  designs  showed 
improved  mixing  and  combustion  efficiency,  but  they  also  significantly  increased 
the  cold  flow  pressure  loss  of  the  flameholder.  The  best  advanced  model  was  ADV5, 
which  used  the  axially  injected  fuel  stream  to  reduce  the  pressure  loss  while 
maintaining  discrete  flame  kernels  at  the  center  of  each  vortex.  Possible  design 
improvements  consist  of  lengthening  the  tail  to  reduce  the  horning  angles  of  the 
ramp,  and  reducing  the  maximum  vertical  extent  of  the  flameholder.  Such 
modifications  will  reduce  the  effectiveness  of  the  vortex,  but  should  make  a  larger 
reduction  to  the  pressure  loss. 

Detailed  combustion  analyses  of  the  triangular  flameholder  were  performed  at 
equivalence  ratios  of  0.61  and  0.55.  Coupling  between  acoustics  and  combustion 
caused  the  flow  to  be  dominated  by  axial  pressure  waves.  The  flow  patterns  were 
irregular  and  characterized  by  the  simultaneous  shedding  of  vortices  from  the  top 
and  bottom  of  the  base  of  the  flameholder.  Simulations  of  a  simpler  flat  plate 
flameholder  at  the  same  equivalence  ratios  showed  the  acoustic-combustion 
interactions  were  reduced  at  the  lower  equivalence  ratio,  but  the  influence  of 
flameholder  geometry  was  not  investigated  further. 

The  isothermal  experimental  tests  were  completed  for  the  baseline  IFF  with  no 
injection.  Regular  alternating  vortices  were  shed  from  the  trailing  edge  with  a 
Strouhal  number  of  0.243,  yielding  a  mean  recirculation  bubble  with  a  length  of  0.9 
H.  The  turbulence  in  the  wake  was  not  isotropic,  with  the  uv  and  w  components 
of  the  Reynolds  stresses  making  the  largest  contribution.  Comparisons  with  the 
numerical  data  showed  that  the  use  of  a  steady-state  symmetric  solution  will 
generate  a  solution  similar  to  a  backward  facing  step  flow,  with  much  lower 
Reynolds  stresses  than  the  bluff  body  flow.  The  use  of  3-D  time-accurate  schemes 
were  necessary  to  properly  capture  the  vortex  shedding  of  the  flowfield.  Such 
algorithms  are  not  practical  for  use  as  design  tools  at  the  present  time. 

A  water-cooled  combustion  test  section  was  designed  and  fabricated.  No  tests  were 
completed  due  to  the  lack  of  availability  of  test  facilities  at  Wright  Laboratory. 
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APPENDIX  A. 

A  NUMERICAL  METHOD  FOR  THE  SECOND-MOMENT 
TURBULENCE  CLOSURE  WITH  APPLICATIONS 
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A.1  INTRODUCTION 


In  the  last  two  decades,  a  large  nximber  of  turbulent  flow  cases  have  been  compiled 
demonstrating  that  the  widely  used  isotropic  eddy-viscosity  closures  is  inaccurate. 
Most  of  these  flow  can  be  classified  as  "complex"  since  extra  significant  rate-of-strain 
components  are  present  due  to  rotation,  streamline  curvature,  separation,  three- 
dimensionality,  etc.  As  a  result,  we  have  witnessed  the  development  of  many 
"improved"  eddy-viscosity  models  to  take  a  particular  extra  strain  rate  into  account 
by  modifying  the  model  constants  or  adding  extra  terms.  Unfortunately,  this 
approach  often  finds  limited  use  due  to  the  lack  of  generality  and  its  ad  hoc  nature. 
A  potentially  more  rewarding  approach  is  to  go  beyond  the  eddy-viscosity 
assumption  and  to  develop  higher  moment  turbulence  closures.  The  second- 
moment  closure  is  the  lowest  level  of  closures  which  has  the  potential  to  retain 
many  important  physical  mechanisms  and  thereby  is  likely  to  account  for  different 
extra-strain  rates  automatically.  It  is  precisely  this  reason  that  the  second-moment 
closure  has  been  the  subject  of  extensive  studies  in  the  last  two  decades.  The 
validation  and  further  refinement  of  newly  developed  second-moment  turbulence 
closures  that  heavily  rely  on  the  use  of  an  accurate  and  robust  computational 
approach.  However,  the  incorporation  of  the  second-moment  closure  into  the 
existing  eddy- viscosity  turbulence  model  based  numerical  approach  is  not  a  trivial 
task.  One  of  the  prominent  issues  is  the  treatment  of  the  Reynolds  stress  terms  in 
the  momentum  equations.  This  treatment  has  a  lot  to  do  with  the  stability  and 
efficiency  of  the  computational  procedure. 

The  major  objective  of  the  present  study  is  to  address  several  numerical  issues 
related  to  the  incorporation  of  the  second-moment  turbulence  closure  into  a 
colocated  grid  approach  on  non-orthogonal  grids.  To  demonstrate  the  robustness  of 
the  numerical  method,  the  Gibson-Launder  model^o  is  incorporated  and  used  to 
calculate  three  turbulent  flows:  backward  facing  step  flow;  faired  diffuser  flow;  and 
confined  swirling  flow. 


A-2 


A.2  GOVERNING  EQUATIONS 


This  study  deals  with  stationary  and  incompressible  turbulent  flows.  The  Reynolds- 
averaged  mass  and  momentum  conservation  equations  can  be  concisely  written  in 
Cartesian  tensor  form  as: 


a{puy) 


=  0 


(1) 


-w- 


(2) 


where  Uj  and  Uj  are  the  jth  components  of  the  mean  and  fluctuating  velocity,  P  is 
the  mean  pressure,  p  and  \i  are  the  fluid  density  and  viscosity.  With  a  second- 
moment  closure,  the  transport  equations  for  the  Reynolds  stresses  are  solved  and 
can  be  written  as: 


(3) 


or  symbolically  as: 


=  +  +  (4) 

where  the  terms  on  the  right  hand  side  (RHS)  represent  the  viscous  diffusion, 
turbulent  diffusion,  production,  redistribution,  and  dissipation  of  iqu}- .  Of  the  five 

terms  on  the  RHS  of  Equation  (3),  Dfj  and  Pij  are  exact  and  require  no  modeling. 

The  remaining  terms  either  involve  higher-order  correlation  of  Ui,  correlations 
with  the  fluctuating  pressure  or  correlations  with  the  velocity  gradient  du/dxj . 

Therefore,  models  are  needed  for  dJ,  and  ejy  to  form  a  closed  set  of  governing 

equations. 

In  this  study,  the  second-moment  closure  follows  closely  the  model  of  Gibson- 
Launder.40  This  selection  is  based  on  the  fact  that  the  Gibson-Launder  model  is  one 
of  the  most  popular  models  and  there  exists  a  large  number  of  validated  cases 
available  for  comparison. 

The  models  for  ufj,  Oij  and  Efy  are  give  as: 
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(5) 

^ij  =  ^ij,l  +  +  ^ijw 

(6) 

%i  =  -Qpf 

(7) 

=  “Q  (P ij  -  -jP 

(8) 

+  Qzy  (^hn.2  “  f  ~  f 

(9) 

£,y  —  5jy 

(10) 

where  nj  is  the  ith  component  of  the  unit  normal  to  a  wall  and 

Jw  -R-  £y^ 

(11) 

with  yw  being  the  distance  normal  to  a  wall. 

Finally,  a  transport  equation  for  the  turbulent  dissipation  rate,  £ 
closure  and  is  modeled  as: 

is  needed  for  the 

^{p^/^)_  a  (.,dE\,  a  frnW 

(12) 

The  model  constants  are  chosen  according  to  Gibson-Launder^o  as: 

{Cp  Cz,  Cg,  C^p  C^,  C^,  k)  =  {1.8, 0.6, 05, 0.3, 0.22, 0.18, 1.44, 1.92,  0.09, 0.41} 
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A.3  NUMERICAL  METHOD 


As  mentioned  in  the  introduction,  the  Reynolds  stress  force,  Bu^j/dxj ,  appears 

explicitly  in  the  momentum  equations  for  the  second-moment  closure.  A  major 
portion  of  this  force  acts  as  a  diffusive  force.  The  numerical  treatment  of  this  force 
needs  special  attention. 

In  the  past,  most  researchers  adopted  the  apparent  anisotropic  eddy  viscosity 
approach.  The  full  Reynolds  stress  force  is  taken  as  a  diffusive  contribution  by 
introducing  an  anisotropic  eddy-viscosity  Pij  as: 

This  way,  the  positive  gij  will  serve  as  an  anisotropic  eddy  viscosity  and  negative  Pij 
will  be  treated  explicitly.  Despite  its  benefit  of  being  able  to  account  for  the 
anisotropic  turbulent  diffusion,  and  ease  of  implementation,  this  approach  may  find 
numerical  problems  associated  with  the  possible  huge  value  of  Pij.  Therefore,  some 
limits  have  to  be  applied  to  Pij  in  the  iterative  procedure  of  a  computational 
method.  Consequently,  the  overall  numerical  procedure  may  be  degraded 
tremendously.  In  this  study,  a  different  approach  is  used.  The  approach  follows  the 
work  of  Obi  et  alM  by  using  a  special  evaluation  of  the  cell  face  Reynolds  stresses  to 
eliminate  the  possible  decoupling  between  velocities  and  Reynolds  stresses  on  a 
colocated  grid.  The  proposed  special  interpolation  technique  resulted  in  an  explicit 
relation  between  the  Reynolds  stresses  and  the  strain  rates.  This  feature  naturally 
leads  to  the  semicoupled  treatment  of  the  Reynolds  stress  force  thereby  promoting 
the  stability  and  robustness  of  the  numerical  method.  The  numerical  method  is 
described  below. 

A.3.1  Discretization 

All  the  governing  equations  except  the  continuity  can  be  expressed  in  Cartesian 
coordinates  as: 


9<{> 


^  d 
3^ 


P  3<t> 
3(f) 


^3K 


+  S, 
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'  3(|)' 


Note  that  Fkk  is  not  a  tensor  and  no  summation  on  T^k  is  implied. 


(14) 
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Equation  (14)  can  be  transformed  into  a  body-fitted  coordinate  (^i,  ^2,  ^3)  as: 


where  is  not  a  tensor  and  is  defined  as: 


+  S* 


and 


d{Xi,X2,X3) 


(15) 


(16) 


(17) 


In  the  above,  e.g.  means  d^j/dxi . 

The  discretization  of  Equation  (15)  is  carried  out  using  a  finite-volume  approach. 
First,  the  solution  domain  is  divided  into  a  finite  number  of  discrete  volumes  or 
cells,  where  all  variables  are  stored  at  their  geometric  centers.  The  equation  is  then 
integrated  over  all  the  cells  by  using  the  Gaussian  Theorem.  A  detailed 
discretization  process  was  described  before  and  is  not  reported  here.  It  is  sufficient  to 
point  out  that  the  final  discretized  equation  for  <|)  at  a  control  volume  P  can  be 
written  in  a  linear  equation  form  as: 

+  ^sw4>sw  +  ^Se4>SE  +  ^NW^NW  Qgs 

+  ^LS^LS  +  ^LN^LN  +  ^H^HS  +  ^HN^HN 
+  ^  M.4>VVL  +  +  ^EL^EL  +  ^EH^EH  +5*^ 

or 

~  ^  ^nb^nb  '^((»  (19) 

In  the  above,  subscripts  W,  E,  S,  N,  L,  H  denote  the  cells  located  West,  East,  South, 
North,  Low,  High  side  of  cell  P,  and  nb  refers  to  all  the  neighboring  cells. 
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In  a  discretized  form,  the  cell  face  stresses  are  needed  to  obtain  the  Reynolds  stress 
force  in  the  momentrun  equation.  To  avoid  the  decoupling  of  the  velocities  and  the 
stresses  and  to  extract  the  necessary  diffusive  contribution  of  the  stress  force,  a 
special  technique  is  needed  and  presented  next. 

First,  the  discretized  Reynolds  stress  transport  equation  is  organized  as  follows: 

_  [  3/j.  dU; 


In  the  above,  V  is  the  cell  volume  and  o^j  is  not  a  tensor  and  no  summation  on  i  or 
j  is  implied.  For  the  Gibson-Launder  model  without  the  inclusion  of  the  wall 
reflection  term,  the  coefficients,  coij  can  be  listed  as: 


Uf  =  (l-§C2jpiJI  = 

y(xiij  =  il-C2)puj 

Note  that  the  above  coefficients  are  always  positive  since  C2  =  0.6  is  used  by  the 
Gibson-Launder  model.  Therefore,  the  cell  center  stress  u^-  can  be  calculated  from 

Equation  (20)  as: 


dUf 

/«"5x7 


(22) 


with 


(23) 


Now,  the  cell  face  Reynolds  stress  can  be  obtained  by  linearly  interpolating  all  terms 
in  Equation  (22)  except  for  the  velocity  gradient,  i.e. 


where  <>  denotes  linear  interpolation  from  cell  centers  to  the  cell  face.  Equation 
(24)  was  used  by  Obi  et  al.  to  calculate  the  cell  face  stresses.  However,  Equation  (24) 
can  be  simplified  by  noting  that: 
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(25) 


Substitution  of  Equation  (25)  into  Equation  (24)  yields: 


(26) 


For  practical  purposes,  it  is  more  convenient  to  use  the  following  interpolation: 


(27) 


It  is  seen  from  Equation  (27)  that  the  cell  face  stress  is  obtained  with  linear  average  of 
the  cell  center  stresses  plus  a  third  derivative  velocity  term.  The  above 
interpolation  formula  is  not  only  much  simpler  to  implement  but  it  also  has  a  clear 
physical  meaning.  The  third-derivative  velocity  term  serves  as  a  high  order 
damping  to  suppress  the  oscillatory  mode.  This  Reynolds  stress  force  in  the 
momentum  equation  has  been  transformed  into  a  new  form  which  gives  an  explicit 
diffusive  contribution. 


In  Body-Fitted  Coordinates,  the  Reynolds  stress  force  can  be  expressed  as: 


Integration  over  a  control  volume  gives: 

=  J  fi^'^  ^j)  (28) 

where  S  denotes  a  central  difference  operator. 

Substitution  of  (27)  into  (29)  yields: 


Ff  =  F,i  + 


(30) 


Pil  = 


idUA  , 


(31) 
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hi 


It  is  seen  by  comparing  Equation  (32)  with  (15)  that  F2i  term  is  exactly  the  diffusion 
term  needed  to  promote  the  stability  of  the  numerical  method  by  defining  for 
Ui-momentum  equation  as: 


(33) 


where  the  superscript,  i,  is  added  to  ^  to  represent  Ui-momentum  equation.  The 
remaining  force  term  fii  will  be  calculated  explicitly. 

A.3.3  Boimdary  Conditions 

Most  boundary  conditions  can  be  imposed  without  any  difficulty  except  at  a  wall  and 
at  a  symmetry  in  a  body-fitted  coordinate  grid.  At  a  solid  wall,  the  wall  function 
approach  of  Launder  and  Spalding  is  used  for  the  momentum  equations.  The 
application  of  the  wall  functions  to  the  Reynolds  stress  equations  is  not  trivial  in  a 
BFC  grid.  In  this  study,  fixed  Reynolds  stress  values  near  a  waU  are  used. 

Naot  et  alA^  derived  the  following  near  wall  stress  relations  based  on  the 
equilibrium  assumptions  and  by  comparing  the  experimental  data. 

^  =  1.95  ;  -p-  =  0-375  ;  *  =  0.675 ;  -  ■?  =  0.2  (34) 

Air  Air  Air  K 
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To  apply  the  boundary  condition  to  a  BFC  grid,  let  us  define  a  local  coordinate 
system  (x,y,z)  at  the  wall  such  that  y  is  normal  to  the  wall;  x  is  along  the  flow 
direction;  and  z  is  normal  to  the  xy  plane.  Then,  the  Reynolds  stresses  defined  on 
the  xyz  system  can  be  obtained  by  Equation  (34).  However,  a  global  cartesian 
coordinate  system  (xi,  X2,  X3)  is  used  in  tiiis  study  and  the  Reynolds  stresses  IT^j  are 

represented  in  (xi,  X2,  X3)  system.  Therefore,  a  transformation  needs  to  be 
established  and  it  can  be  expressed  as: 

UjTj  =  aijOiyU^  +  Ct2f^2p'  +  «3i«3 ^  ;)  ^  (35) 


where  aij  is  the  directional  cosine  representing  the  projection  of  the  ith  unit  base  of 
(x,y,z)  system  to  the  jth  base  of  the  (xi,  X2,  X3)  system.  With  Equation  (35),  the  wall 
boundary  conditions  can  be  easily  applied  to  the  stress  equations  in  a  general 
coordinate  system. 
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On  a  symmetry  surface,  a  similar  transformation  is  applied.  The  Reynolds  stresses 
on  a  symmetry  surface  can  be  obtained  as: 


with 

^lm= 

where  is  the  stress  at  the  cell  center. 


(36) 


(37) 
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A.4  APPLICATIONS 


Three  turbulent  flows  are  calculated  to  demonstrate  the  numerical  method.  The 
first  is  the  backward  facing  step  flow  investigated  by  Eaton  and  Johnston.  43  This  case 
is  used  to  make  sure  the  BFC  formulation  recovers  to  the  cartesian  grids  as  required. 
The  second  case  is  a  turbulent  flow  over  a  faired  annular  diffuser  measured  by 
Stevens  and  Fry.  44  This  case  is  chosen  to  demonstrate  the  BFC  capability.  The  third 
case  is  the  confined  swirling  flow  measured  by  Nejad  et  alA5  and  is  used  to 
demonstrate  the  solution  of  all  six  Reynolds  stress  equations.  The  second-moment 
calculation  results  will  be  compared  to  the  standard  k-e  model  to  show  the  strength 
of  the  second-moinerit  closure. 

A.4.1  Backward-Facing  Step  Flow 

The  computational  domain  is  shown  in  Figure  A-1.  The  inlet  U- velocity  and  the 
turbulent  kinetic  energy  are  obtained  from  the  experimental  data.  The  V-velocity 
component  and  the  Reynolds  shear  stresses  are  set  to  zero  and  the  normal  Reynolds 
stresses  are  obtained  assuming  isotropy.  The  dissipation  rate  at  the  inlet  is  estimated 
from  Ki-5/(0.05H).  Computations  are  performed  on  a  nonvmiform  grid  consisting  of 
140x70  cells  with  40  streamwise  cells  located  between  x/H  =  (-1.5, 1.0)  and  another  60 
between  x/H  =  (1.0,  12.0).  This  grid  should  be  sufficient  to  generate  a  grid- 
independent  solution  based  on  the  previous  results.  The  measure  of  the 
computational  time  indicated  that  the  second-moment  closure  took  about  4  to  5 
more  total  CPU  time  than  the  k-e  model. 


Figure  A-1.  The  Computational  Domain  of  the  Backward-Facing  Step 

The  mean  streamwise  velocity  profiles  are  presented  in  Figure  A-2  comparing  the 
second-moment  closure  and  the  k-e  model  with  the  experimental  data.  Similar  to 
the  finding  of  Obi  et  al.  where  a  different  backward  facing  step  was  calculated,  the 
second-moment  closure  is  able  to  predict  the  experimentally  observed  monotonic 
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increase  in  the  reverse  velocity  away  from  the  bottom  wall  (see  x/H  =  1.0  in  Figure 
A-2),  whereas  the  k-e  model  behaves  in  an  opposite  trend.  This  different  behavior 
is  due  to  the  fact  that  the  second-moment  closure  predicts  the  existence  of  a 
secondary  recirculation  at  the  corner  while  the  k-e  model  fails  to  predict  it.  This  can 
be  seen  clearly  in  Figure  A-3  where  the  computed  streamlines  are  displayed.  The 
predicted  size  of  the  secondary  recirculation  is  about  0.93H  which  is  quite  consistent 
with  the  experimental  data  of  l.OH.  Overall,  it  is  seen  that  the  second-moment 
closure  gives  better  mean  velocity  agreement  with  the  experimental  data.  Despite 
quite  different  predictions  of  the  secondary  recirculation,  the  calculated  sizes  of  the 
primary  recirculation  are  similar  between  the  two  models.  The  second-moment 
model  predicts  the  reattachment  length  of  7.06H  while  the  k-e  model  is  6.91H.  The 
experimentally  determined  reattachment  length  is  8.0H.  Therefore,  both  models 
underpredict  the  reattachment  length  by  12-14%. 


O  :Data,Eaton&Jolmstoiil980 

-  :  Sec  ond-  M  om  ent  M  odel 

-  ■  :K-£ Model 
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Figure  A-3.  The  Calculated  Streamlines  -Top:  the  second-moment  closure 

Bottom:  the  k-e  model 


Figures  A-4  and  A-5  present  the  turbulent  normal  stress  and  the  shear  stress  W  . 

Both  stresses  are  predicted  well  by  the  second-moment  closure  at  x/H  =  1.0  and  2.0, 
but  greater  discrepancies  show  up  further  downstream.  As  expected,  the  normal 
stress  is  severely  underpredicted  by  the  k-e  model  due  to  its  isotropic  eddy-viscosity 
assumption.  However,  the  shear  stress  is  calculated  unexpectedly  well  by  the  k-e 
model  as  seen  in  Figure  A-5. 

A.4.2  An  Axisymmetric  Faired  Diffuser  Flow 

The  turbulent  flow  through  an  axisymmetric  faired  diffuser  as  shown  in  Figure  A-6 
is  calculated  to  demonstrate  the  calculation  method  on  a  body-fitted-coordinate  grid. 
The  diffuser  consists  of  an  entry  annular  pipe,  an  inlet  bend,  a  straight  diffuser,  an 
outlet  bend,  and  an  exit  annular  pipe.  The  experimental  test  was  carried  out  by 
Stevens  and  Fry  and  a  numerical  study  was  reported  by  Jones  &  Manners.  The  inlet 
of  the  computational  domain  is  located  at  0.0762  m  upstream  of  the  inlet  bend. 
(Note  that  the  width  of  the  entry  annular  pipe  is  0.0254  m.)  All  main  variables  at 
the  inlet  are  obtained  from  the  fully-developed  annulus  solutions  to  minimize  the 
inlet  condition  uncertainty.  This  is  possible  since  a  long  entry  length  (about  50.0 
hydraulic  diameter  of  the  entry  annulus)  was  used  in  the  experiment  and  the  data 
suggested  that  the  fully  developed  conditions  have  been  established. 


o 


Data,  ii:aton&  Johnston  1980 
Seconnd-M  omentM  odei 


—  —  —  :K-£  Model 


Figure  A-4.  The  Comparison  of  the  Turbuleiit  Normal  Stress  Component 
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-  :  Second-Moment  Model 

—  —  —  ;S-£  Model 


Figure  A-5.  The  Comparison  of  the  Turbulent  Shear  Stress  Component 


Figure  A-6.  The  Geometry  of  the  Axisymmetric  Faired  Diffuser 
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The  exit  of  the  computational  domain  is  placed  0.508  m  downstream  of  the  outlet 
bend  and  all  the  variables  are  extrapolated  except  the  pressure.  A  fixed  pressure 
condition  is  applied  at  the  exit.  The  long  length  of  the  exit  annulus  is  used  to 
minimize  the  influence  of  the  Newman  type  exit  boundary  conditions. 

Both  the  second-moment  turbulence  model  of  Gibson-Launder  and  the  standard  k-e 
model  are  used  for  the  calculation.  The  computational  mesh  is  chosen  as  100x50  in 
the  streamwise  and  cross-stream  directions,  respectively,  based  on  the  grid- 
independence  study  of  Jones  and  Manners.  The  second-moment  model  took  about 
twice  the  total  CPU  time  of  the  k-e  model  calculation.  Both  computations  will  reach 
the  convergence  criterion  in  600  to  800  iterations. 

Figures  A-7  and  A-8  present  the  radial  distribution  of  the  mean  axial  velocity  and 
Reynolds  shear  stress  at  several  streamwise  locations,  respectively.  The  station 
numbers  correspond  to  the  selected  measurement  stations  as  indicated  in  Figure 
A-6.  Ti  is  the  local  nondimensional  distance  coordinate  approximately  normal  to 
the  inner  and  outer  annular  walls  with  t]  =  0  at  the  outer  wall.  As  can  be  seen,  from 


the  outer  wall  and  underpredicted  at  the  inner  wall  by  the  k-e  model.  This  trend  of 
the  results  by  the  k-e  model  is  not  surprising  since  the  streamline  is  strongly  curved 
for  the  flow.  The  flow  streamline  curvature  tends  to  invalidate  the  standard  k-e 
model.  Therefore,  the  calculated  result  prove  that  the  second-moment  closure  has 
naturally  taken  the  extra-strain  effect,  induced  by  streamline  curvature,  into 
account.  The  most  remarkable  result  is  the  comparison  of  the  streamwise  velocity 
at  station  8  in  Figure  A-7.  The  experimental  data  suggest  that  the  maximum 
velocity  location  is  switched  from  the  outer  wall  to  the  inner  wall  after  the  fluid  is 
coming  out  of  the  outlet  bend.  This  unusual  feature  is  well  captured  by  the  second- 
moment  closure.  However,  the  k-e  model  will  keep  the  maximum  velocity 
location  near  the  outer  wall,  and  totally  missed  this  experimentally  observed 
feature. 


Figure  A-7,  the  second-moment  closure  prediction  of  the  mean  streamwise  velocity 
In  the  straight  diffuser  section,  the  mean  streamwise  velocity  is  overpredicted  near 
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Figure  A-7.  The  Comparison  of  the  Mean  Streamwise  Velocity 
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Figure  A-8.  The  Comparison  of  the  Turbulent  Shear  Stress  Component 
A.4.3  Confined  Swirling  Flow 

The  confined  swirling  flow  calculated  in  this  study  was  experimentally  studied  by 
Nejad  et  al.^s  and  is  shown  schematically  in  Figure  A-9.  It  is  seen  that  the  dump 
combustor  consists  of  two  sections:  the  inlet  pipe  and  the  combustion  chamber. 
The  inlet  pipe  has  an  iimer  diameter  of  101.6  mm  and  a  cylindrical  teflon  swirler  is 
located  50.8  mm  upstream  of  the  dump  plane.  The  combustor  chamber  is  1850  mm 
in  length  and  has  a  152.4  mm  I.D.  Air  is  the  working  fluid  and  the  inlet  centerline 
velocity  is  measured  360  mm  upstream  of  the  swirler  housing  and  is  kept  at  19.2 
m/s.  This  corresponds  to  a  Reynolds  number  of  125,000  based  on  the  combustor 
inlet  diameter.  The  computation  corresponds  to  the  case  with  a  swirl  number  of  0.5. 
The  experimental  data  showed  that  this  swirl  number  is  just  strong  enough  to  lead 
to  flow  reversal  near  the  chamber  centerline  after  the  dump  plane. 
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Figure  A-9.  The  Geometry  of  the  Swirling  Flow  Studied 

For  all  the  computations,  the  inlet  is  located  0.38  step  height  downstream  of  the 
dump  plane.  This  is  the  closest  location  to  the  swirler  where  detailed  experimental 
data  are  available  for  the  mean  velocity  components  and  the  Reynolds  stresses 
except  e.  It  is  known  that  swirling  flow  computations  are  very  sensitive  to  the  inlet 
conditions.  Therefore,  complete  availability  of  the  experimental  data  at  the  inlet 
makes  this  set  of  data  valuable  to  assess  turbulence  model  performance.  The  only 
uncertainty  at  the  inlet,  therefore,  is  the  specification  of  the  turbulent  dissipation 
rate  e  and  it  will  be  discussed  shortly. 

The  outlet  boimdary  is  placed  72.8  step  heights  downstream  of  the  dump  plane. 
This  outlet  is  far  enough  from  the  inlet  to  minimize  the  influence  of  the  exit 
boundary  conditions.  At  the  outlet,  all  main  variables  are  extrapolated. 

At  the  solid  walls,  the  standard  wall  function  approach  is  adopted.  Finally,  at  the 
centerline  symmetry,  the  derivatives  of  all  variables  are  set  to  zero  except  ®  and 
m  •,m  =  im  =  Q  at  symmetry  instead. 

The  turbulent  dissipation  rate  must  be  specified  by  the  modeler.  One  of  the  most 
used  approaches  is  based  on  the  turbulence-energy  equilibrium  and  e  is  estimated  as 

y1-5 

e  =  ^  (38) 

where  R  is  the  radius  of  the  chamber  and  a  is  a  constant  to  be  determined.  In  the 
RSTM  calculation  of  strongly  swirling  flows  by  Hogg  and  Leschziner46,  a  of  0.36  was 
used  based  on  numerical  optimization.  In  another  k-e  computation  of  the  Nejad  et 
al.45  swirling  flow,  a=0.3  was  used  by  Favaloro  et  al.47  since  it  gave  the  best 
compromise  for  several  swirl  numbers  calculated.  In  this  study,  an  independent 
calculation  was  carried  out  to  investigate  the  sensitivity  of  the  solution  to  a.  Two  a- 
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values,  0.36  and  1.8,  were  used  to  calculate  the  swirling  flow  using  the  RST  model. 
The  calculated  radial  distributions  of  the  mean  axial  velocity,  swirl  velocity,  and  the 
turbulent  kinetic  energy  are  shown  in  Figure  A-10  at  two  axial  locations.  It  is  seen 
that  the  calculated  results  are  quite  sensitive  to  the  inlet  e.  It  is  suggested  that  the 
inlet  e  should  be  chosen  such  that  the  predicted  resultant  turbulent  kinetic  energy 
matches  the  experimental  data  as  close  to  the  inlet  as  possible.  The  above  criterion 
is  based  on  two  observations:  (1)  k  directly  responds  to  the  inlet  e  and  is  heavily 
influenced  by  e.  Besides,  k  is  readily  measurable  and  is  available  from  most 
turbulent  flow  data  sets.  (2)  Theoretically,  6  should  be  chosen  such  that  all  the 
measured  data  satisfy  the  governing  equations  at  the  inlet.  This  is  obviously 
impractical.  Therefore,  the  nearest  data  to  the  inlet  could  be  used  to  determine  the 
proper  inlet  8.  In  this  study,  the  above  criterion  results  in  the  choice  of  a  of  0.36  as 
used  by  Hogg  and  Leschziner.46  It  should  be  pointed  out  in  passing  that  the  same 
inlet  £  also  resulted  in  close  agreement  with  the  experimental  data  at  x/h=1.0  for  the 
k-e  model. 

In  this  study,  effort  has  been  spent  to  use  the  grid  which  provides  essentially  a  grid- 
independent  solution.  Two  highly  nonuniform  grids  are  selected:  a  70x35  grid 
(designated  as  coarse  grid)  and  a  120x60  grid  (designated  as  fine  grid).  Most  of  the 
axial  grid  is  clustered  within  the  central  and  behind-step  separated  zones.  For  the 
fine  grid,  30  cells  are  used  between  the  inlet  and  x/h  of  4,  and  60  cells  are  non- 
uniformly  distributed  between  x/h  of  4  and  x/h  of  26.  The  RSTM  predicted  radial 
variations  of  the  mean  axial  velocity,  swirl  velocity,  and  the  turbulence  kinetic 
energy  are  displayed  in  Figure  A-11  at  two  axial  locations.  It  is  seen  that  increasing 
grid  from  75x35  to  120x60  produced  only  marginal  differences.  Therefore,  it  is 
believed  that  the  solutions  obtained  using  the  120x60  grid  are  essentially  grid 
independent.  It  is  worthy  to  mention  that  a  48x48  grid  was  selected  by  Hogg  and 
Leschziner‘16  and  a  100x66  grid  was  used  by  Jones  and  Pascau^^  based  on  their  grid- 
independence  study. 

The  Reynolds  stress  transport  model,  the  standard  k-e  model,  and  the  RNG-based  k- 
e  model  are  used  to  calculate  the  swirling  flow  of  Nejad  et  al.4S  All  the  results 
presented  below  are  obtained  using  the  fine  grid  (120x60).  The  results  are  considered 
to  be  converged  only  if  the  total  absolute  sum  of  the  residuals  for  each  governing 
equation,  normalized  by  the  inlet  flux,  is  reduced  four  orders  of  magnitude.  Further 
decrease  in  the  residuals  did  not  change  the  final  solution  more  than  0.1  percent. 
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O  C  :  ExTDerimentaiData 

- - :  Prediction  with  a=0,36 

_ :  Prediction  with  a=:1.3 


O  G  :  Experimental  Data 

- —  :  Prediction  with  a=sOG6 

- -  :  Prediction  with  a=L3 


(a)  Mean  Axial  Velocity 


(b)  Mean  Swirl  Velocity 


(c)  Turbulent  Kinetic  Energy 


Figure  A-10.  Sensitivity  of  the  Computed  Solution  to  the  Inlet  e 
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(b)  Mean  Swirl  Velocity 
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(c)  Turbulent  Kinetic  Energy 


Figure  A-11.  Sensitivity  of  the  Computed  Solution  to  the  Grid 
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The  calculated  streamline  patterns  of  the  three  turbulence  models  are  displayed  in 
Figure  A-12  and  the  predicted  size  of  the  central  separation  zone  (Lc)  and  the  comer 
reattachment  length  (Lt)  are  compared  to  the  experimental  data  in  Table  A-1.  It  is 
seen  from  Figure  A-12  and  Table  A-1  that  both  the  standard  k-e  model  and  the 
RNG-based  k-e  model  failed  to  predict  the  existence  of  the  central  reversed  zone.  On 
the  other  hand,  the  RST  model  was  capable  of  predicting  the  occurrence  of  the 
central  reversed  zone  although  there  was  an  over-prediction  of  the  size.  It  is 
believed  that  the  failure  of  the  two-equation  type  models  is  resulted  from  the 
isotropic  eddy-viscosity  assumption.  This  isotropy  leads  to  a  too-high  turbulent 
viscosity  for  the  swirling  flow  and  consequently  suppresses  the  central  reversed 
zone. 


Table  A-1.  Comparison  of  the  Central  Separation  Size  (Lc)  and  the  Corner 
Reattachment  Length  (Lt) 


Computational  Results 

Experimental  Data  j 

Lc/H 

Lt/H 

Lc/H 

Lt/H 

RST  Model 

5.96 

3.15 

4.4 

32 

Standard  k-c 
Model 

0.84 

3.53 

4.4 

32 

RNG-Based 
k-e  Model 

0.89 

4.36 

4.4 

3.2 

To  substantiate  the  above  argument,  the  flow  was  also  calculated  using  the  k-e 
model  with  much  lower  inlet  [It-  This  could  be  achieved  by  increasing  the  inlet 
turbulence  dissipation  rate  e.  It  was  found  that  the  central  reversed  flow  would 
appear  if  inlet  e  was  several  times  larger.  As  expected,  however,  the  calculated 
turbulence  quantities  such  as  the  turbulent  kinetic  energy  deviated  too  far  away 
from  the  experimental  data  near  the  inlet  section.  This  points  to  the  inadequacy  of 
the  use  of  high  inlet  8.  These  results  suggested  that  the  Reynolds  stress  transport 
model,  which  implicitly  uses  the  anisotropic  turbulent  viscosity,  seems  to  have  the 
basic  mechanisms  capable  of  more  accurately  predicting  the  central  separation  of 
swirling  flows.  The  overprediction  of  the  separation  size  by  the  RST  model  may  be 
attributable  to  the  use  of  Gibson-Launder  pressure-strain  model.  Computations 
with  more  advanced  pressure-stain  models,  such  as  those  by  Shih  et  al.49  and 
Speziale  et.  al.50^  could  clarify  the  above  argument  but  are  not  within  the  scope  of 
this  study. 
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(a)  The  RST  Model 


(c)  The  RNG-Based  k-e  Model 


Figure  A-12.  The  Calculated  Streamline  Plots 
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The  prediction  of  the  corner  recirculation  zone  is  less  sensitive  to  the  inlet 
condition.  Overall,  the  RST  model  provides  the  best  agreement  to  the  experimental 
data  for  the  comer  reattachment  length.  It  is  seen  that  both  eddy-viscosity  models 
over  predicted  the  reattachment  length,  a  behavior  quite  contrary  to  the  case  of  a 
sudden  expansion  flow  without  swirl.  Experimentally,  it  was  found  that  the 
reattachment  decreased  from  approximately  8  step  heights  for  the  nonswirl  flow  to 
about  3.2  step  heights  for  S  =  0.5  case.  This  decrease  in  the  comer  recirculation  is 
caused  by  the  rapid  expansion  of  the  flow  induced  by  the  action  of  the  centrifugal 
forces.  Tlierefore,  the  above  results  suggested  that  the  eddy-viscosity  model  fails  to 
take  the  action  of  the  swirl  generated  centrifugal  force  into  account  properly.  This 
may  be  directly  linked  to  the  failure  of  the  model  to  predict  the  existence  of  the 
central  separation  zone.  It  should  point  out,  in  passing,  that  the  RNG-based  k-e 
model  consistently  predicted  a  larger  comer  recirculation  zone  than  the  standard  k- 
e  model  no  matter  if  the  flow  has  swirl  or  not.  If  no  swirl  is  present,  the  RNG 
model  tends  to  have  a  good  prediction  of  the  reattachment  length  (see  Speziale  and 
Thangam29)  since  the  standard  k-e  model  underpredicted  the  reattachment  length. 
As  shown  in  this  study,  however,  the  RNG  model  can  predict  worse  results  for  the 
sudden  expansion  flow  with  swirl. 

Comparisons  of  the  calculated  and  measured  radial  variation  of  the  mean  axial 
velocity  and  the  swirl  velocity  are  shown  in  Figure  A-13  and  A-14  respectively,  at 
several  axial  locations.  Again,  the  RST  model  has  the  best  overall  agreement  with 
the  experimental  data.  The  results  between  the  standard  and  RNG-based  k-e  models 
are  essentially  similar.  For  the  mean  axial  velocity  component,  the  major  difference 
between  the  results  of  the  RST  model  and  that  of  the  eddy-viscosity  model  are 
confined  to  the  central  separation  zone.  This  is  a  direct  consequence  of  the  failure  of 
the  eddy-viscosity  model  to  predict  the  reversed  flow  at  the  center.  The  most 
noteworthy  difference  between  the  results  of  different  models  is  the  prediction  of 
the  mean  swirl  velocity  far  downstream  of  the  dump  plane  (x/h  >  5.0).  As  is  clear  in 
Figure  A-14,  the  two  eddy-viscosity  models  predicted  that  the  swirling  flow  will 
eventually  approach  to  the  solid-body-rotation  type.  This  solid-body-rotation  flow 
will  remain  until  the  end  of  the  calculation  domain.  (This  is  the  reason  x/h  of  18.0 
station  is  added  to  Figure  A-14).  However,  the  experimental  data  suggested  that  the 
flow  retained  the  strength  of  its  vortex  core  all  the  way  to  the  exit  and  did  not 
evolve  to  solid-body-rotation  flow.  This  interesting  feature  was  captured  quite  well 
by  the  RST  model.  The  above  results  indicated  that  the  RST  model  is  able  to 
account  for  the  complicated  interaction  between  swirl  and  turbulence  field  and  to 
predict  the  combined  presence  of  free  and  forced  vortex  in  the  swirling  flow. 
However,  the  eddy-viscosity  model  tends  to  lose  the  strength  of  the  vortex  cone  and 
expand  to  a  solid-body  type  rotational  motion. 
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Figure  A-13.  Radial  Variation  of  the  Mean  Axial  Velocity  at  Four  Axial  Locations 
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Figure  A-14.  Radial  Variation  of  the  Mean  Swirl  Velocity  at  Five  Axial  Locations 

Finally,  Figure  A-15  shows  that  comparison  of  the  calculated  and  measured 
turbulence  kinetic  energy.  As  discussed  in  the  previous  section,  the  turbulence 
kinetic  energy  was  predicted  quite  well  by  all  models  near  the  inlet  by  choosing  a  = 
0.36  for  the  inlet  dissipation  rate,  (x/h  =  1.0  for  the  present  case).  But  further 
downstream,  the  eddy-viscosity  models  predicted  much  lower  turbulence  kinetic 
energy  levels.  This  underprediction  results  from  the  improper  modeling  of  the 
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swirl  and  turbulence  interaction  by  the  eddy-viscosity  model  and  is  probably 
responsible  for  the  failure  of  the  model  to  preserve  the  strength  of  the  vortex  core 
and  to  predict  the  reversed  flow  at  the  centerline.  Again,  the  RST  model  provided 
the  best  overall  prediction.  However,  large  discrepancies  could  be  noticed, 
particularly  around  x/h  of  5.0  and  near  the  centerline  at  downstream  locations.  The 
underprediction  of  k  around  x/h  of  5.0  is  probably  responsible  for  the  overprediction 
of  the  central  separation  zone.  With  limited  information  available,  no  definite 
conclusion  can  be  drawn  with  regard  to  this  discrepancy.  One  major  source  of  error 
could  come  from  the  pressure-strain  model. 


r/h 


Figure  A-15.  Radial  Variation  of  the  Turbulent  Kinetic  Energy  at  Four  Axial 
Locations 
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A.5  CONCXUSIONS 


The  second- moment  turbulence  closure  is  implemented  into  a  finite-volume 
pressure  based  numerical  method  on  colocated  and  nonorthogonal  grids.  Special 
interpolation  scheme  is  developed  to  remove  the  decoupling  between  momentum 
and  Reynolds  stress  equations  and  to  provide  the  strong  diffusive  contribution  of 
the  Reynolds  stress  force.  Three  two-dimensional  turbulent  flows  are  selected  to 
demonstrate  the  developed  numerical  method  and  the  performance  of  the  second- 
moment  closure.  The  comparison  between  the  second-moment  closure  and  the  k-e 
model  with  the  experimental  data  finds  that  the  second-moment  closure  does  offer 
advantages  over  the  k-e  model  for  complex  turbulent  flows.  For  some  cases,  the 
difference  between  the  model  performances  could  be  remarkable.  Due  to  the 
increased  number  of  partial  differential  equations  to  be  solved,  the  second-moment 
turbulence  will  consume  more  computer  resources.  For  the  three  flows  calculated 
in  this  study  the  increase  in  CPU  time  for  the  second-moment  closure  is  between 
two  to  six  times  that  of  the  k-e  model. 
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APPENDIX  B. 

LASER  DOPPLER  VELOOMETRY 
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B.l  LDV  THEORY 


Laser  Doppler  Velocimetry  (LDV)  has  been  used  for  fluid  flow  measurement  since 
the  mid  1960's.  Durst,  et  al.si  describe  the  theory  and  application  of  LDV  at  length. 
LDV  is  a  nonintrusive  instrument  used  to  measure  turbulent  velocity  fields 
without  disturbing  the  flow,  and  is  used  where  conventional  techniques  (i.e.,  hot 
wire)  are  not  practical.  The  only  added  requirement  for  LDV  data  acquisition  is  the 
presence  of  seed  material  in  the  flow  to  illvuninate  the  laser  beams. 

The  basic  operational  principle  of  the  LDV  system  stems  from  the  creation  of  a  probe 
volume  by  crossing  two  polarized  laser  beams.  At  the  crossing,  the  two  beams 
constructively  and  destructively  interfere  to  create  fringes  of  light  as  shown  in 
Figure  B-1.  As  the  seed  particles  in  the  flow  pass  through  the  probe  volume,  the 
reflected  light  from  the  fringes  generate  an  amplitude  modulated  sinusoidal  signal, 
called  a  Doppler  burst.  The  scattered  light  is  collected  by  a  lens  and  focused  onto  a 
photomultiplier,  which  transmits  an  electrical  signal  to  the  signal  processor.  The 
LDV  signal  processors  amplify  and  filter  the  signals  from  the  photomultiplier, 
validate  the  Doppler  frequency,  and  compute  the  Doppler  period  which  is  the 
reciprocal  of  the  Doppler  frequency.  The  frequency  of  the  scattered  light  as  the 
particle  crosses  bright  and  dark  fringes  is  related  to  tire  particle  velocity  component 
perpendicular  to  the  fringe  direction  by  the  following  relation: 


FOCUSING  LENS 
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Figure  B-1.  Typical  Dual  Beam  Laser  Doppler  Anemometer 
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Notice  that  the  velocity  is  a  linear  function  of  the  Doppler  frequency,  since  all  other 
terms  are  constant  for  a  given  LDV  configuration. 

One  shortcoming  of  LDV  is  that  a  particle  will  produce  the  same  frequency 
regardless  of  the  direction  of  travel.  Therefore  a  particle  traveling  at  40  m/s  will 
produce  the  same  frequency  as  one  traveling  at  -40  m/s  since  there  is  not  a  negative 
frequency.  This  is  shown  graphically  in  Figure  B-2,  and  is  known  as  aliasing.  If 
negative  velocities  are  expected,  the  problem  can  be  solved  by  imposing  a  frequency 
shift  on  the  probe  volume.  The  frequency  shift  causes  the  fringes  in  the  probe 
volume  to  move  at  a  constant  rate.  A  particle  sitting  at  rest  within  the  probe 
volume  will  then  register  a  Doppler  frequency  equal  to  the  frequency  shift  (typically 
40  MHz).  This  has  the  effect  of  shifting  the  curve  up  the  ordinate  in  Figure  B-2  by  40 
MHz.  The  curve  again  reflects  about  f=0,  and  aliasing  will  again  occur  for  particles 
traveling  less  than  the  base  velocity  (Uq).  The  maximum  velocity  is  limited  by  the 
accuracy  of  the  photomultiplier  and  the  signal  processor,  and  is  limited  by  f^ax- 

The  frequency  range  can  also  be  shifted  in  the  opposite  direction,  to  increase  the 
peak  velocity  that  the  system  can  measure.  This  is  often  done  to  measure 
supersonic  flows.  In  short,  there  is  a  range  of  velocity  acquisition  that  falls  between 
f=0  and  f=fimx,  and  that  range  can  be  moved  along  the  velocity  scale  through  the  use 
of  frequency  shifting. 

The  frequency  shifting  can  accomplished  by  several  methods,  one  of  which  is  to  pass 
one  of  the  laser  beams  through  an  acoustic  oscillator  known  as  a  Bragg  cell.  The 
Bragg  cell  creates  a  diffraction  grating  using  the  air  density.  As  the  beam  passes 
through  the  grating,  the  mean  is  split  into  an  infinite  number  of  beams,  each  shifted 
in  frequency  by  an  integer  value  of  the  diffraction  grating.  The  first  order  shifted 
beam  is  separated  from  the  rest  and  is  focused  to  cross  with  the  original  beam  and 
form  the  probe  volume. 

Two  and  three  dimensional  LDV  systems  are  created  by  focusing  additional  pairs  of 
laser  beams  of  different  colors  at  the  same  location,  but  at  different  orientations. 
Each  pair  of  beams  must  have  its  own  photomultiplier,  with  appropriate  filters  so 
that  only  the  correct  color  light  is  registered. 
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Figiire  B-2.  Frequency  Shifting  Allows  the  Measurement  of  Negative  Velocities 
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B.2  STATISTICAL  ANALYSES 


LDV  systems  measure  many  samples  of  individual  velocity  realizations.  The  most 
meaningful  form  of  the  data  comes  from  standard  statistical  analysis  of  the 
measured  data.  Since  there  are  inherent  errors  involved  in  statistical  analysis  based 
on  sample  population,  the  statistics  are  considered  estimates,  and  not  actual 
measurements.  A  list  of  the  statistical  equations  used  in  the  current  study  are  as 
follows: 


Mean: 

or 

RMS: 

^rms=y 

or 

n 

V  n 

SIGMA: 


/  n  \  f  n 


n 


or 


where  n  is  the  number  of  samples  in  the  data  set. 


Turbulence  intensity  (%):  100  Ua/Uref 


Second  Moments: 
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Triple  Products: 
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APPENDIX  C 

FORTRAN  LISTING  FOR  LDV  DATA  REDUCTION 
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The  following  FORTRAN  program  was  written  to  allow  the  proper  reduction  of  the 
LDV  data.  The  program  allows  for  the  correction  in  any  angle  discrepancies  that 
may  be  present  in  a  nearly  orthogonal  laser  system.  Point-by-point  normalization  of 
the  data  is  also  possible.  The  program  takes  up  to  two  input  files  in  addition  to  the 
raw  data  files  produced  by  the  data  acquisition  program.  These  files  are: 

family.PRN  -  Header  file  contairung  reference  conditions  for  each  point. 

The  data  must  be  ordered  to  match  the  numerical  order  of 
the  raw  data  file. 

master  -  Optional  file  that  contains  all  input  data  to  the  program 
such  as  the  family  name,  number  of  files  to  reduce,  and 
correction  angles.  This  allows  many  files  of  different 
families  to  be  reduced  in  a  batch  mode.  If  the  file  does  not 
exist,  the  user  is  prompted  for  each  input. 

A  sample  header  file  is  presented  in  Table  C-1  and  is  self  explanatory.  The  only 
value  actually  used  by  the  program  is  the  reference  velocity.  The  coordinates  are  not 
cross-matched  to  the  raw  data  file.  Therefore,  the  entries  in  the  header  file  must 
correspond  consecutively  with  data  file  00,  01,  02,  etc.  Even  if  the  data  files  do  not 
start  with  00,  the  header  file  must  start  with  an  entry  for  00. 

Two  sets  of  angles  are  input  to  the  program.  The  first  set  is  the  correction  angles  to 
produce  an  orthogonal  system  between  the  blue  and  green  beams  (Figure  C-1).  The 
included  angle  between  the  beams  should  always  be  90°  +  6b  -  0g.  Therefore,  there  is 
an  infinite  number  of  correct  combinations  of  9b  and  0g  for  a  given  configuration. 

The  second  set  of  angles  is  the  rotation  and  tilt  of  the  LDV  system.  Rotation  (5)  is 
defined  as  angle  between  the  green/blue  beams  and  the  reference  axes  (Figure  C-2). 
The  tilt  (e)  is  defined  as  the  rotation  of  the  entire  LDV  system  about  the  x-axis 
(Figure  C-3). 

For  a  typical  case,  9g  =  0,  0b  =  included  angle  -  90, 

5  =  angle  between  the  green  beams  and  the  tunnel  x-axis 
e  =  0 
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Figure  C-1.  Definition  of  Correction  Angles  (0g,  0b) 


C-3 


tunnel  vertical  coordinate 


violet  transmitter  head 


green/blue  transmitter  head 


tunnel  spanwise  coordinate 


Figure  C-3.  Definition  of  Tilt  Angle  (e) 


c  program  name  is:round5.for 
c  date:  5f3H94 

c  program  does  data  reduction  for  as  many  files  in  master 
c  program  prints  raw  velocity  data  with  header  (.y**) 
c  it  calculates  statistics  from  raw  velocity  data  (.sta) 
c  correlation  coefficients,  3rd  &  4th  moments, 
c  skewness  &  flatness  coefficients  (.mis) 
c  it  makes  proplot  files  which  are  quantities  normalized 
c  (.pll):  velocities,  sigmas  &  turbulent  stresses 
c  (.pl2):  triple  products 

c  (.pl3):  normal  stresses  such  as  uu,vv,ww,uuu,vvv,www 
c  execute  program  from  tsiXdata 
c  header  is  in  tsiXdata 

implicit  double  precision  (a-h,o-z) 
characterise  master,angl,answer 
characterise  filer,family,familyn 
integer  first,last,npts,nsig 
logical  inter 
integeri2  proc2,proc3 

realiS  mom3u(5e),mom3v(5e),mom3w(5e),mom4u(5e),mom4v(5e), 
&mom4w(5e),vel(12eee,3),tbd(12eee) 
integer  head,raw,yfile,sta,mis,pll  ,pl2,pl3 
common/stat/Uiu,Viu,Wiu,Ui2,Vi2,Wi2,UiVi, 
&UiWi,ViWi,UiViWi,UiUiWi,ViViWi,UiViVi,UiUiVi,UiWiWi,ViWiWi, 
&Ui3,Vi3,Wi3,Ui4,Vi4,Wi4,tt,ncount 
common  /  unit/  y  file, is  ta,imis,ipll  ,ipl2,ipl3,iy  file 
conunon/  angle/  anglu,anglv,anglw 
dimension  hl(5e),h2(5e),h3(5e),h4(5e),h5(5e),h6(50),h7(5e) 
dimension  ubar(5e),vbar(50),wbar(5e),uv(5e),uw(5e),vw(5e) 
dimension  usig(5e),vsig(5e),wsig(50) 
dimension  uvw(5e),uuw(50),uuv(5e),uww(5e),uvv(5e), 

&vvw(5e),vww(5e) 

dimension  cocouv(5e),cocouw(5e),cocovw(5e) 
dimension  flatu(5e),flatv(5e),flatw(5e) 
dimension  skewu(5e),skewv(5e),skeww(5e) 
common/stat2/h4,h5,h6,h7,ubar,vbar,wbar,usig,vsig,wsig, 
&uv,uw,vw,uvw,uuw,uuv,uww,uvv,vvw,vww,cocouv,cocouw, 
&cocovw,flatu,flatv,flatw,skewu,skewv,skeww, 
&mom3u,mom3v,mom3w,mom4u,mom4v,mom4w 
mas=e5 
head=ie 
raw  =  15 
yfile=ll 
sta=21 
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inis=41 

pll=51 

pl2=52 

pl3=53 

22  print*,'please  enter  master  name:(c/r  for  interactive)  ' 
read(*,24)master 
if  (master.eq.'  ')then 
inter=.true. 
else 

inter =. false. 

call  opener(mas,master/  '/old  ’derr) 
endif 

if(ierr.ne.O)goto  22 

c  OPEN(99,nLE='JUNK\JUNK3.DAr) 
do  296  nm=l,100 

if(inter)print^'please  enter  family  name  (c/ r  to  end):  ' 
read(mas,24)  family 
ifffamily.eq.'  ')stop 

if(inter)print*/please  enter  first  file  #:  ’ 
read(mas,*)first 

if(inter)prmt*/please  enter  last  file  #  (.a.  for  all):  ' 
read(mas,24)answer 
if  (answer  .eq. '  a'  )then 
last=100000 
else 

read(answer/)last 

endif 

if(inter)PRINT*, 

&  'please  enter  new  family  name  (c/r  for  old  name):  ' 
read  (mas,24)  familyn 

24  format  (3a) 

if  (familyn.eq.'  ')  familyn =family 

if(inter)print*/please  enter  #  of  standard  of  deviations:  ' 
read(mas,'*)nsig 
c 

anglu=0. 

anglv=0. 

anglw=0. 

if(inter)print*, 

&:'please  enter  correction  angles:(c/r  for:0,0,0)  ' 
read(mas,24)angl 
if  (angleq.'  ’)goto  25 
read(angl,*)anglu,anglv,anglw 

25  pi=acos(-l.) 
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anglu=anglu*pi/ 180. 
anglv=(anglv+90)*pi/ 180. 
anglw=anglw*pi/ 180. 
nfiles=last-first+l 
c  '.pm'  is  the  header  file 

call  opener(head,family,'.pm'/old  ',ierr) 

if(ieiT.ne.0)stop 

readChead,"^) 

readChead,"*) 

readQiead,’^) 

readOhead,’^) 

do  29  i=l,&st 

read(head/) 

29  continue 

c  this  section  reads  &  prints  all  in  one  row  in  "e-format" 
if(inter)print*, 

&'please  enter  (y)  if  you  want  .sta  files:(c/r  for  none)’ 
read(mas,24)answer 
if(answer.ne.'y')  goto  77 
ista=l 

call  opener(sta,familyn,'.sta','unknown'4err) 
if  (ierr.ne.O)  stop 

write(sta,42)  ’X’/Y'/Z’/Ts’/Ps'/Dp'/Uref’/U/Uref’, 
&'V/Uref,'W/Uref’/Usig/Uref/Vsig/Uref'/Wsig/Uref', 
&’UV/Uref^2’;UW/Uref^2*,'VW/Uref^2','lJU/Uref''2’;W/UreP2', 
&'WW/Uref^2’,'UVW/UreP'3'/UUV/UreP3','UUW/Uref^3';VVU/Uref^3', 
&'VVW/Uref^3';WWU/Uref^3’;WWV/Uref^3';UUU/Uref^3*, 
&'WV/Uref'^3’;WWW/Uref^3';UV-CORR’;UW-CORR’,’VW-CORR', 
&'U-3rd'/V-3rd';W-3rd';U-Skew','V-Skew';W-Skew', 
&'U-4th’,'V-4th','W-4th';U-Flat’;V-Flat';W-Flat' 

42  format(lX,A3,4X,A3,3X,A3,2X,A3,4X,A3,2X,A3,7X,A5,8X,a6, 
&9x,a6,8x,a6,5x,a9,2(5X,A9),2(5x,a9),2(5X,A9),2(5x,a9), 
&5x,al0,6(4X,A10),4x,al0,4x,al0,4x,al0, 

&6X,A7/6x,A7,7x,a7,9x,a5,9x,a5,10X,A5, 

&7X,A6,8X,A6,8X,A63(9x,a5), 

&9X,A6,7x,a6,8x,a6) 

write(sta,53)  'in'/in'/in';R'/psi','in-H20','m/s', 
&'COEf';coef';coef;mom';mom';mom’;coef;coef', 
&'COef';mom';mom';mom’,'COef';coef';coef' 

53  fonnat(lX,A3,4X,A3,3X,A3,3X,A2,3X,A4,lX,A6,6X,A3, 
&100x400x,118x,a4,9x,a4,10x,a4,12x,a3,llx,a3,12x,a3,9x,a4, 
&2(10x,a4)^(10x,a4),llx,a4,9x,a4,10x,a4) 

77  if(inter)print*, 

(Sc'please  enter  (y)  if  you  want  .mis  files:(c/ r  for  none)' 
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read(inas,24)answer 
if(answer.ne.'y')  goto  96 
imis=l 

c  this  section  prints  corr.  coeff.,  moments,  &  misc. 
call  opener(mis,familyn/.mis'/unknown',ierr) 
if  (ierr.ne.O)  stop 
write(mis/) 
write(inis,’^) 
write(mis,*) 
write(mis,*) 

write(mis,89)  ’X'/Y'/Z’/Ts’/Ps'/Dp'/Uref/CoUV, 

&'CoUW’;CoVW’;3rd’;3rd’;3rd’;SkewU';SkewV’, 

&'SkewW’;4th';4thV4thVFlatU’;FlatV';FlatW' 

89  fonnat(lX,A3,4X,A3,3X,A3,3X,A2,2X,A4,lX,A4,3X,A5, 
&4x,a4,4x,a4,4x,a4,6X,A3,6X,A3,6X,A3,5x,a5,4x,a5,4x, 
&a5,6X,A3,7X,A3,7X,A3,4X,A5,3X,A53X,A5) 
write(mis,95)  'inVin’,'in'/R','psiVin-H20', 
i&'m/s'/Coef’/Coef'/Coef/MomU’/MomV'/MomW’/Coef’, 
&:'Coef’/Coef;MomU’;MomV’;MomW’;Coef';Coef’;Coef' 

95  format(lX,A3,4X,A3,3X,A3,3X,A2,2X,A4,lX,A6,2X,A3, 
&5x,a4,4x,a4,4x,a4,6X,A4,5X,A4,5X,A4,5x,a4,5x,a4,5x,a4, 
&6X,A4,6X,A4,6X,A43x,a4,4x,a4,4x,a4) 

96  if(inter)print*, 

(Sc’please  enter  (y)  if  you  want  .pll  files:(c/r  for  none)’ 
read(mas,24)answer 
iffanswer.ne.'y')  goto  99 
ipll=l 

c  this  section  prints  the  proplot  (.pll)  files 

call  opener(pll,familyn,’.pll’/unknown',ierr) 
if  (ierr.ne.O)  stop 

write(pll,98)  'X'/Y’/Z’/U/Uref'/V/Uref/W/Uref, 
&’Usig/Uref';Vsig/Uref’;Wsig/Uref’;UV/UreP2', 
&’UW/Uref'^2';VW/UreP2’ 

98  format(lX,A3,4X,A3,3X,A3,3(6X,a6),3(4X,A9),3(4X,A9)) 
c  this  section  prints  the  proplot  (.pl2)  files 

99  if(inter)print*, 

&'please  enter  (y)  if  you  want  .pl2  files:(c/r  for  none)' 
read(mas,24)answer 
if(answer.ne.’y')  goto  110 
ipl2=l 

call  opener(pl2,familyn,'.pl2',’unknown',ierr) 
if  (ierr.ne.O)  stop 

write(pl2,100)  'X'/Y’/Z'/UVW/Uref^O’/UUV/Uref'^O', 
&’UUW/Uref^3';VVU/Uref^3';VVW/Uref^3';WWU/Uref'^3’, 
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&'WWV/Uref^3' 

100  fonnat(lX,A3,4X,A3,3X,A3,7(4X,al0)) 
c  this  section  prints  the  proplot  (.pl3)  files 

110  if(inter)print*, 

&’please  enter  (y)  if  you  want  .pl3  files:(c/r  for  none)' 

read(mas,24)answer 

iffanswer.ne.’y')  goto  108 

ipl3=l 

call  opener(pl3,familyn,'.pl3’,'unknown',ierr) 
if  (ierr.ne.O)  stop 

write(pl3433)  'X’/Y'/Z’/UU/Uref'^Z/VV/Vref^T, 
&'WW/Wref'^2';UUU/Uref^3’;WV/Uref'^3'/WWW/Uref^3’ 
133  format(lX,A3,4X,A3,3X,A3,6(4X,al  0)) 
c  this  section  prints  the  (.y**)  files 

108  if(inter)print*, 

&'please  enter  (y)  if  you  want  .y**  files:(c/r  for  none)' 
read(mas,24)answer 
if  (answer. ne.'y')  goto  111 
iyfile=l 

111  do  256  n=l,nfiles 

call  rdfhdr(raw,family,firstn,npts,vel,tbd,proc2,proc3,del,eps) 
read(head,%end=257)Hl(n),H2(n),H3(n),H4(n),H5(n),H6(n),H7(n) 
if(iyfile.ne.l)goto  125 
writeffiler,!  09)first+n-l 
if(first+n-l.ge.l0)write(filer,112)first+n-l 
call  opener(yfile,fainil)m,filer,'unknown',ierr) 
if  (ierr.ne.0)  stop 

109  fonnat('.y0',il) 

112  formate. y',i2) 

write(yfile,118)  'X(in)’/Y(in)’;Z(in)';Ts(R)’;Ps(psi)', 
&'Dp(in-H20)';Uref(m/s)' 

118  format(5X,A7,5X,A7,4X,A7,4X,A7,4X,A8,4X,All,4X,All) 
write(yfile,120)hl(n),h2(n),h3(n),h4(n),h5(n),h6(n),h7(n) 

120  format(6X,F5.2,7X,F5.2,6X,F5.2,7X,F5.1,4X,F8.4,8X,F4.2,9X,F8.5) 
writeCyfile,*) 

write(yfile424)  'Ni’/Ui(m/s)';Vi(m/s)’,’Wi(m/s)';TBD(s)', 
&;’TT(s)' 

124  format(3X,A5,8X,A8,6X,A8,6X,A840X,A8,10X,A8) 

125  stdu=lE28 
stdv=lE28 
stdw=lE28 
ubar(n)=0 
vbar(n)=0 
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wbar(n)=0 

ksig=2 

if(nsig.eq.O)  then 

ksig=l 

endif 

iysave=iyfile 
do  236  isig=l,ksig 
iyfile=iysave+isig-ksig 

call  statsuin(vel,tbd,stdu,stdv,stdw,proc2,proc3,npts,n) 

print*,ncount 

Ubar(n)=Uiu/ncount 

Vbar(n)= Viu/ ncount 

Wbar(n) = Wiu  /  ncoimt 

UV(n)=UiVi/  ncount-Ubar(n)*Vbar(n) 

UW  (n)=UiWi  /  ncount-Ubar(n)*Wbar(n) 

VW(n)=ViWi /ncount- Vbar(n)*Wbar(n) 
Usig(n)=sqrt(Ui2/ncount-Ubar(n)*Ubar(n)) 
Vsig(n)=sqrt(abs(Vi2/ncount-Vbar(n)*Vbar(n)))+(l-proc2) 
Wsig(n)=sqrt(abs(Wi2 /ncount- Wbar(n)*Wbar(n)))+(l-proc3) 
STDU=NSIG=^USIG(n) 

STDV=NSIG*VSIG(n) 

STDW=NSIG’^WSIG(n) 

236  continue 

if  (ista+ipl2.ge.  1  )then 

UVW(n)=UiViWi/ncount+2'^Ubar(n)*Vbar(n)*Wbar(n)-Vbar(n)* 
&UiWi/ncount-Ubar(n)*ViWi/ncount-Wbar(n)*UiVi/ncount 
UUV(n)=UiUiVi/ncount+2*Ubar(n)’^Ubar(n)*Vbar(n)-Vbar(n)* 
&Ui2/ ncount-2*Ubar(n)*UiVi/ ncount 
UUW(n)=UiUiWi/ncount+2*Ubar(n)*Ubar(n)*Wbar(n)-Wbar(n)'^ 
&Ui2/ ncount- 2*Ubar(n)*UiWi/ ncount 
UVV(n)=UiViVi/ncount+2’*Ubar(n)*Vbar(n)*Vbar(n)-Ubar(n)* 
&Vi2/ncount-2*Vbar(n)*UiVi/ncount 
VVW(n)=ViViWi/ncount+2’^Vbar(n)*Vbar(n)’^Wbar(n)-Wbar(n)* 
&Vi2/ncount-2’^Vbar(n)*ViWi/ncount 
UWW(n)=UiWiWi/ncount+2*Ubar(n)*Wbar(n)*Wbar(n)-Ubar(n)* 
&Wi2/ncount-2*Wbar(n)*UiWi/ncount 
VWW(n)=ViWiWi/ncount+2*Vbar(n)*Wbar(n)*Wbar(n)-Vbar(n)* 
&Wi2/ncount-2*Wbar(n)*ViWi/ncount 
endif 

if(ista+imis.ge.l)then 

MOM3U(n)=Ui3/ ncount-3*Ui2/ ncounb*Ubar(n)+2*Ubar(n)**3 
MOM3V(n)=Vi3/ncount-3’^Vi2/ncount*Vbar(n)+2*Vbar(n)**3 
MOM3W(n)=Wi3/ncount-3*Wi2/ncount*Wbar(n)+2’^Wbar(n)**’3 
SKEWU(n)=MOM3U(n)  /(2=^Usig(n)*=^3) 
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SKEWV(n)=MOM3V(n)/  (2*Vsig(n)**3) 

SKEWW(n)=MOIvI3  W(n)  /  (2*Wsig(n)**3) 

MOM4U(n)=Ui4/ncount-4*Ui3/ncount*Ubar(n)+6*Ubar(n)**2* 

&Ui2/ncount-3*lJbar(n)**4 

MOM4  V  (n) = Vi4  /  ncoiint-4*Vi3  /  ncount*  Vbar(n) +6  *Vbar(n)**2* 

&  Vi2  /  ncount-3*Vbar(n)*’*’4 

MOM4W(n)=Wi4/ncount-4*Wi3/ncount*Wbar(n)+6*Wbar(n)**2* 

&Wi2/ncount-3*Wbar(n)*M 

FLATU(n)=MOM4U(n)/USIG(n)*M 

FLATV(n)=MOM4V(n)/VSIG(n)’^’^4 

FLATW(n)=MOM4W(n)/WSIG(n)*’»4 

COCOUV(n)=UV(n)/(Usig(n)*Vsig(n)) 

COCOUW(n)=UW(n)  /  (Usig(n)*Wsig(n)) 
COCOVW(n)=VW(n)/(Vsig(n)*Wsig(n)) 
endif 

256  continue 

257  nfiles=n-l 
if(hl(l).ne.hl(2))  then 

call  sortGil,h2,h3,nfiles) 
else  if(h2(l).ne.h2(2))  then 
call  sort(h2,hl,h3,nfiles) 
else 

call  sortQi3,hl,h2,nfiles) 
endif 

c  this  section  prints  all 
do  260  n=l,nfiles 
if  (ista.ne.l)goto  243 

write  (sta,241)  hl(n),h2(n),h3(n),h4(n),h5(n),h6(n),h7(n), 
&ubar(n)/h7(n),vbar(n)/h7(n),wbar(n)/h7(nXusig{n)/h7(n), 
&vsig(n)/h7(n),wsig(n)/h7(n),uv(n)/h7(n)*%uw(n)/h7(n)**2, 
&vw(n)/h7(n)*-^,usig(n)’^*2/h7(n)**2,vsig(n)**2/h7(n)**2, 
&wsig(n)**2/h7(n)**2,uvw(n)/h7(n)**3,uuv(n)/h7(n)**3, 

&UU  w(n)  /  h7(n)**3,uvv(n)  /  h7(n)**3,vvw(n)  /h7(n)**3,uww(n)  /h7(n)**3, 
&vww(n)/h7(n)*=^3,usig(n)=»*3/h7(n)'^*3, 

&vsig(n)’^*3/h7(n)*’^,wsig(n)**3/h7(n)**3,cocouv(n),cocouw(n), 

&cocovw(n),mom3u(n),mom3v(n),inom3w(n),skewu(n),skewv(n),skeww(n), 

&moin4u(n),mom4v(n),mom4w(n),flatu(n),flatv(n),flatw(n) 

241  fonnat(lx,f5.2,lx,f5.2,lx,f5.2,lx,f5.1,lx,f5.2,lx,f4.2, 
&4(lx,el3.6),3(lx,el3.6),3(lx,el3.6),7(lx,el3.6), 

&21(lx,el3.6)) 

c  this  section  prints  moments  and  misceilaneous(.mis) 

243  if(imis.ne.l)goto  244 
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write  (mis,249)  hl(n),h2(n),h3(n),h4(n),h5(n),h6(n),h7(n), 
&cocouv(n),cocouw(n),cocovw(n), 

&:mom3u(n),moin3v(n),mom3w(n),skewu(n),skewv(n),skeww(n),mom4u(n), 

&inom4v(n),mom4w(n),flatu(n),flatv(n),flatw(n) 

249  format(lx,f5.2,lx,f5.2,lx,f5.2,2x,i3,lx,f5.2,lx,f4.2, 
&4(lx,f7.4)^(lx,f8.2),3(lx,f8.3),3(lx,f9.2),3(lx,f7.3)) 

c  this  section  prints  proplot  (.pll) 

244  if(ipll.ne.l)goto  245 

write  (pll,250)  hl(n),h2(n),h3(n),ubar(n)/h7(n),vbar(n)/h7(n), 
&wbar(n)/h7(n),usig(n)/h7(n), 

&vsig(n)  /  h7(n),  wsig(n)  /  h7(n),uv(n)  /  (h7  (n)*h7  (n)),uw(n)  /  (h7(n)’^ 
&h7(n)),vw(n)/  (h7(n)*h7(n)) 

250  format(lx,f5.2,lx,f5.2,lx,f5.2,3(lx,ell.4)^(lx,el2.5), 

&3(lx,el2.5)) 

c  this  section  prints  proplot  (.pl2) 

245  if(ipl2.ne.l)  goto  246 

write  (pl2,251)  hl(n),h2(n),h3(n),uvw(n)/h7(n)**3,uuv(n)/h7(n)**3, 

&uuw(n)  /  h7  (n)  **3, 

&uvv(n)/h7(n)**3,vvw(n)/h7(n)’^*3,uww(n)/h7(n)**3,vww(n)/h7(n)**3 

251  fonnat(lx,f5.2,lx,G.24x,f5.2,7(lx,el3.6)) 
c  this  section  prints  proplot  (.pl3) 

246  if(ipl3.ne.l)goto  260 

write  (pl3,252)  hl(n),h2(n),h3(n),usig(n)’^*2/h7(n)**2, 

«Sc  vsig(n)  **2  /h.7(n)**2, 

&wsig(n)**2/h7(n)**2,usig(n)**3/h7(n)**3,vsig(n)**3/h7(n)’*’*3, 

&wsig(n)*’^/h7(n)**3 

252  fornaat(lx,f5.2,lx,f5.2,lx,f5.2,6(lx,el3.6)) 

260  continue 

296  continue 
stop 
end 

subroutine  sort(x,y^^pts) 

implicit  double  precision  (a-h,o-2) 

dimension  x(*),y(*),z(*) 

dimension  h4(50),h5(50),h6(50),h7(50) 

dimension  ubar(50),vbar(50),wbar(50) 

dimension  usig(50),vsig(50),wsig(50),uv(50),uw(50),vw(50) 

dimension  uvw(50),uuw(50),uuv(50),uww(50),uvv(50), 

&:w  w(50),  vww(5  0) 

dimension  cocouv(50),cocouw(50),cocovw(50) 
dimension  flatu(50),flatv(50),flatw(50) 
dimension  skewu(50),skewv(50),skeww(50) 

real*8  mom3u(50),mom3v(50),mom3w(50),mom4u(50),mom4v(50). 
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&moin4w(50) 

cominon/stat2/h4,h5,h6,h7,ubar,vbar,wbar,usig,vsig,wsig, 
&uv,uw,vw,uvw,uuw,uuv,uww,uvv,vvw,vww,cocouv,cocouw, 
&cocovw,flatu,flatv,flatw,skewu,skewv,skeww, 
&mom3u,inom3v,mom3w,moin4u,mom4v,mom4w 
do  10  i=l,npts-l 
do  10  j=i+l,npts 
if  (x(j).ltx(i))then 
temp=x(i) 
x(i)=x(j) 
x(j)=temp 
temp=y(i) 
y(i)=y(j) 
y(j)=temp 
temp=z(i) 
z(i)=z(j) 
z(j)=temp 
temp=h4(i) 
h4(i)=h4(j) 
h4(j)=temp 
temp=h5(i) 
h5(i)=h5(j) 
h5(j)=temp 
temp=h6(i) 
h6(i)=h6(j) 
h6(j)=teinp 
temp=h7(i) 
h7(i)=h7(j) 
h7(j)=temp 
temp=ubar(i) 
ubar(i)=ubar(j) 
ubar(j)=temp 
temp=vbar(i) 
vbar(i)=vbar(j) 
vbar(j)=temp 
temp=wbar(i) 
wbar(i)=wbar(j) 
wbar(j)=temp 
temp=usig(i) 
usig(i)=usig(j) 
usig(j)=temp 
temp=vsig(i) 
vsig(i)=vsig(j) 
vsig(j)=teinp 


temp=wsig(i) 

wsig(i)=wsig(j) 

wsig(j)=temp 

teinp=uv(i) 

uv(i)=uv(j) 

uv(j)=temp 

teinp=uw(i) 

uw(i)=uw(j) 

uw(j)=temp 

teinp=vw(i) 

vw(i)=vw(j) 

vw(j)=temp 

teinp=uvw(i) 

uvw(i)=uvw(j) 

uvw(j)=temp 

teinp=uuv(i) 

uuv(i)=uuv(j) 

uuv(j)=temp 

temp=uuw(i) 

uuw(i)=uuw(j) 

uuw(j)=temp 

temp=uw(i) 

uw(i)=uvv(j) 

uw(j)=temp 

temp=vvw(i) 

ww(i)=vvw(j) 

vvw(j)=temp 

temp=uww(i) 

uww(i)=uww(j) 

uww(j)=temp 

temp=vww(i) 

vww(i)=vww(j) 

vww(j)=temp 

temp=cocouv(i) 

cocouv(i)=cocouv(j) 

cocouv(j)=temp 

temp=cocouw(i) 

cocouw(i)=cocouw(j) 

cocouw(j)=temp 

temp=cocovw(i) 

cocovw(i)=cocovw(j) 

cocovw(j)=temp 

teinp=mom3u(i) 

mom3u(i)=mom3u(j) 


inoin3u(j)=temp 

teinp=mom3v(i) 

moin3v(i)=mom3v(j) 

inom3v(j)=temp 

teinp=mom3w(i) 

mom3w(i)=inom3w(j) 

inom3w(j)=temp 

temp=skewu(i) 

skewu(i)=skewu(j) 

skewu(j)=temp 

temp=skewv(i) 

skewv(i)=skewv(j) 

skewv(j)=teinp 

temp=skeww(i) 

skeww(i)=skeww(j) 

skeww(j)=temp 

teinp=mom4u(i) 

mom4u(i)=mom4u(j) 

mom4u(j)=temp 

temp=mom4v(i) 

mom4v(i)=mom4v(j) 

mom4v(j)=temp 

teinp=mom4w(i) 

mom4w(i)=mom4w(j) 

mom4w(j)=temp 

teinp=flatu(i) 

flatu(i)=flatu(j) 

flatu(j)=temp 

temp=flatv(i) 

flatv(i)=flatv(j) 

flatv(j)=temp 

teinp=flatw(i) 

flatw(i)=flatw(j) 

flatw(j)=temp 

endif 

10  continue 
return 
end 

c  subroutine  rdfhdr 

c  parameters:  unit:  in  -  unit  number  to  read  file  data 
c  description:  the  following  reads  the  file  header  for  a  raw 
c  or  statistics  file 
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subroutine  rdfhdr(unit,fainily,first,nfile,totpts,vel,tbd, 

&  proc2,proc3,del,eps) 

implicit  integer*2  (a-z) 
character*20  drparr(153) 
character*80  filer,family 
character*4  ext 

real=^fshift(3),df(3),vel(12000,3),tbd(12000),fd,del,eps, 
&anglu,anglv,anglw,pi,sinu,sinv,cosu,cosv,cosw,sinuv,sinvu, 
&cosdel,coseps,sindel,sineps,ux,uy,uz,u,v,w 
dimension  hword(4),buffer(1000) 
common/ angle/ anglu,anglv,anglw 
integer*4  temp  l,temp2,hword4(2),tottim 
data  hword/  Z00FF',Z'0300',Z'4000',Z’0FFF7,tottim/0/ 
data  hword4  /  Z’EFFFFFFF’,  Z'OOOOFFFF'/ 
write(ext,  1 09)first+nfile-l 
if(first+nfile-l.ge.l0)write(ext,n2)first+nfile-l 
109  format('.rO',il) 

112  format('.r’,i2) 

call  opener(unit,family,ext,'old  'derr) 
if  (ierr.ne.0)retum 
l=lngth(family) 
write(filer,l  0)family  (1  :l),ext 
10  format(2a) 

writeC*,*) 'processing  file  ’,filer 
do  100  i=l^ 
read(unit,399)drparr(i) 

100  continue 

c  read  header 

do  150  i=6,91 
read(unit,399)  drparr(i) 

150  continue 

c  read  the  dxmuny 

c  information 

do  200  i=92,152 
read(unit,399)  drparr(i) 

200  continue 

read(unit,499)  drparr(153) 
c 

close(unit) 

read(drparr(7),  599)  totadd 
read(drparr(8)  ^99)  tbdflag 
read(drparr(10),599)  nkpts 
read(drparr(24),699)  df(l) 
read(drparr(25),699)  df(2) 
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read(drparr(26),699)  df(3) 
read(drparr(28),699)  fshift(l) 
read(drparr(29),699)  fshift(2) 
read(drparr(30),699)  fshift(3) 
read(drparr(74),799)  ttiine,proc2,proc3 
read(drparr(143),699)del 
read(drparr(144),699)eps 
totword=totadd*(2+ttiine) 
proc2=0 
proc3=0 

if(totadd.eq-2)then 

proc2=l 

proc3=0 

else  if  (totadd.eq.3)then 
proc2=l 
proc3=l 
end  if 

if  (tbdflag  .ge.  1)  totword  =  totword+2 
open(unit,file=filer,access='direct',recl=840) 
nrec=nkpts*1024*totword/  420 
nlast=nkpts*l 024*totword-nrec‘^420 
pi=acos(-l.) 
sinu=sin(anglu) 
smv=sin(anglv) 
cosu=cos(anglu) 
cosv=cos(anglv) 
cosw=cos(anglw) 
sinuv=sin(anglu-anglv) 
sinvu=sin(anglv-anglu) 
del=del*pi/180. 
eps=eps’^i/180. 
cosdel=cos(del) 
coseps=cos(eps) 
sindel=sin(del) 
sineps=sin(eps) 
noff=0 
totpts=0 

do  300  lrec=5,5+nrec 
do  280  ioff=l,noff 

280  buffer(ioff)=buffer(420+noff2-noff+ioff) 
noff2=noff 

if  (lrec.eq.5+nrec)noff2=nlast+noff-420 
read(iimt,rec=lrec)(buffer(i),i=noff+l,420+noff2) 
ipts=(420+noff2)  /  totword 
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noff={420+noff)-ipts*totword 
do  300  jpts=l,ipts 
totpts=totpts+l 
do  275  iadd=0,totadd-l 
iword=(jpts-l)*totword+l+iadd*(2+ttiine) 
ncount=iand(buff er(iword)  ,hword(l ) ) 
if  (ncount.  eq.O )  then 
totpts=totpts-l 
go  to  301 
endif 

nadd=ishft(iand(biiffer(iword),hword(2)),-8)+l 
mant=iand(buffer(iword+l),hword(4)) 
exp=ishft(buffer(iword+l  ),-12) 
Fd=ncovint/(mant’^dble(2.0)**(exp-3)*le-9)*le-6 
vel(totpts,nadd)=(Fd-Fshift(nadd))’*Df(nadd) 

275  continue 

if(tbdflag.gt.O)  then 
tempi  =buffer(iword+2+ttime) 
temp2=buffer(iword+3+ttime) 
templ=ishft(templ,16)+iand(temp2,hword4(2)) 
tbd(totpts)=(templ-tottim) 
tottim=templ 
endif 

300  continue 

301  continue 
tbd(l)=0. 

do  310  i=l,totpts 

Ux=(vel(i,2)*sinu-vel(i,l  )*sin  v)  /  sinuv 
Uy=(vel(i,2)*cosu-vel(i,l  )’^cos  v)  /sinvu 
U2=(vel(i,3)*cosw) 

U=Ux*cosdel-Uy*sindel 

V=Ux*coseps*sindel+Uy'*coseps*cosdel-Uz*sineps 
W=Ux*sineps*sindel+Uy*sineps*cosdei+Uz*coseps 
vel(i,l)=U 
vel(i,2)=V 
vel(i,3)=W 
310  continue 
399  format(a20) 

499  for[nat(al4) 

599  format(i2) 

699  format(g20.0) 

799  fonnat(9x,il,6x,2il) 
return 
end 
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subroutine  opener(imit,root,ext,stat,ierr) 

characterise  root,fname 

character*?  stat 

character*4  ext 

integer  unit 

l=lngth(root) 

write(fname,l  0)root(l  :l),ext 
10  format(2a) 

open  (unit,file=fname,status=stat,iostat=ierr) 
if(ierr.eq.O)  return 

write(i,i)'error  opening  file:’,fname 
return 
end 

fxmetion  Ingth(string) 
character*(i)  string 
do  3  lngth=len(string),l,-l 
if(string(lngth:lngth).ne.'  Oretum 
3  continue 
return 
end 

subroutine  statsum(vel,tbd,stdu,stdv,stdw,proc2,proc3,npts,n) 
integer  yfile 

implicit  double  precision  (a-h,o-z) 
dimension  rset(23),vel(12000,3),tbd(12000) 
common/stat/Uiu,Viu,Wiu,Ui2,Vi2,Wi2,UiVi, 
&UiWi,ViWi,UiViWi,UiUiWi,ViViWi,UiViVi,UiUiVi,UiWiWi,ViWiWi, 
&Ui3,Vi3,Wi3,Ui4,Vi4,Wi4,tt,ncount 
common  /unit/yfile,ista,imis,ipll/ipl2,ipl3,iyfile 
equivalence  (rset(l),Uiu) 
integer*2  proc2,proc3 

dimension  h4(50),h5(50),h6(50),h7(50) 
dimension  ubar(50),vbar(50),wbar(50) 
dimension  usig(50),vsig(50),wsig(50),uv(50),uw(50),vw(50) 
dimension  uvw(50),uuw(50),uuv(50),uww(50),uvv(50), 

&ww(50),vww(50) 

dimension  cocouv(50),cocouw(50),cocovw(50) 
dimension  flatu(50),flatv(50),flatw(50) 
dimension  skewu(50),skewv(50),skeww(50) 
real*8  mom3u(50),mom3v(50),mom3w(50),mom4u(50),mom4v(50), 
&mom4w(50) 
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common/stat2/h4,h5,h6,h7,ubar,vbar,wbar,usig,vsig,wsig, 
&uv,uw,vw,uvw,uuw,uuv,uww,uvv,vvw,vww,cocouv,cocouw, 
&cocovw,flatu,flatv,flatw,skewu,skewv,skeww, 
&mom3u,mom3v,mom3w,mom4u,moin4v,mom4w 
do  15  i=l,23 
rset(i)=0 
15  continue 
ncount=0 
z=0.0e-16 
do  20  i=l,npts 
f3=vel(i,l) 

f5=vel(i,2)+(l-proc2) 

f7=vel(i,3)+(l-proc3) 

f8=tbd(i) 

if(i.eq.l)f8=0 

dfc=3.4399 

dfo=3.4399 

f7c=f7*dfc/dfo 

tt=f8+tt 

if(abs(ubar(n)-f3).gt.stdu)f3=z 
if(abs(vbar(n)-f5).gt.stdv)f5=z 
if(abs(wbar(n)-f7c).gt.stdw)f7c=z 
if  (iyfile.eq.l )  write(yfile,*)i,f3,f5,f7c^8  /  le6,tt/ 1  e6 
if((f3.eq.z).or.(f5.eq.z).or.(f7c.eq.z))go  to  20 
if(abs(f3).lt.le-15)f3=0.0 
if(abs(f5).lt.le-15)f5=0.0 
if(abs(f7C).ltle-15)f7C=0.0 
ncount=ncount  +1 
U=f3 

V=f5-(l-proc2) 

W=f7c-(l-proc3) 

Uiu=Uiu+U 

Viu=Viu+V 

Wiu=Wiu+W 

Ui2=Ui2+(U**2) 

Vi2=Vi2+(V*-^) 

Wi2=Wi2+(W*-^) 

UiVi=UiVi+(U*V) 

UiWi=UiWi+(U*W) 

ViWi=ViWi+(V’»W) 

UiViWi=UiViWi+(U*V*W) 

UiUiWi=UiUiWi+(U’^U*W) 

ViViWi=ViViWi+(V*V*W) 

UiViVi=UiViVi+(U’^V*V) 
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UiUiVi=UiUiVi+(U*U*V) 

UiWiWi=UiWiWi+(U*W*W) 

ViWiWi=ViWiWi+(V’*W*W) 

Ui3=Ui3+(U*=^3) 

Vi3=Vi3+(V’**3) 

Wi3=Wi3+(W’^*3) 

Ui4=Ui4+(U**4) 

Vi4=Vi4+(V*M) 

Wi4=Wi4+(W**4) 

20  continue 
return 
end 

block  data 
integer  yfile 

common  /imit/yfile,ista,imis,ipll,ipl2,ipl3,iytile 
data  ista,imis,ipll4pI2,ipl3,iytile  /  6’*0  / 
end 
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